PN expressions

𝓔(pnsystem)
binding_energy(pnsystem)

Compute the gravitational binding energy of a compact binary.

Note that this may not be as useful as its derivative, 𝓔′, which is used as part of the right-hand side for orbital evolutions.

The nonspinning orbital binding energy is known through 4pN. The expressions through 3.5pN here come from Eq. (233) of Blanchet (2014).

The 4pN term from Eq. (5.2d) of Jaranowski and Schäfer is known exactly, now that the $ν$-linear piece is given as Eq. (32) of Bini and Damour (2013a). The remaining terms are not known exactly, but Bini and Damour (2013b) have derived some terms, though there is incomplete information, which are noted as the constants in this code.

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).

The spin-squared terms (by which we mean both spin-spin and spin-orbit squared terms) in the energy are known to 3pN order, and given in Eq. (3.33) of Bohé et al. (2015).

The spin-cubed terms are known to 3.5pN order, and come from Eq. (6.17) of Marsat (2014).

The tidal-coupling terms come in to the binding 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 $χ₁ = m₁/m$ in their notation; in particular, $χ$ is not a spin parameter. Also note that $λ̂ = λ₂ v^{10}/(m₁+m₂)^5$, and we need to add the coupling terms again with $1 ↔ 2$. Finally, note the normalization difference, where a different overall factor is used, leading to a sign difference.

𝓔′(pnsystem)
binding_energy_deriv(pnsystem)

Compute the derivative with respect to $v$ of the binding energy of a compact binary.

This is computed automatically (via FastDifferentiation) from 𝓔; see that function for details of the PN formulas.

𝓕(pnsystem)
gw_energy_flux(pnsystem)

Compute the gravitational-wave energy flux to infinity

The nonspinning flux terms are complete to 4.5pN order, and are given in Eq. (6.11) of Blanchet et al. (2023).

The spin-orbit terms in the flux are now known to 4.0pN. These terms come from Eq. (4.9) of Marsat et al. (2013)

The spin-squared terms (by which we mean both spin-spin and spin-orbit squared terms) in the flux are known to 3pN order, and given most conveniently in Eq. (4.14) of Bohé et al. (2015).

The spin-cubed terms are known to 3.5pN order, and come from Eq. (6.19) of Marsat (2014).

Beyond 4.5pN, terms are only known in the extreme-mass-ratio limit. These terms are given in Appendix A of Fujita (2012). He computed them up to 22pN. That seems like overkill, so we'll just go up to 6pN.

For systems with matter, the tidal-heating terms come in at relative 5pN order, and are known partially at 6pN order. These terms come from Eq. (3.6) of Vines et al. (2011). Note their unusual convention for mass ratios, where $χ₁ = m₁/m$ in their notation; in particular, $χ$ is not a spin parameter. Also note that $λ̂ = λ₂ v^{10}/(m₁+m₂)^5$, and we need to add the coupling terms again with $1 ↔ 2$. Finally, note the normalization difference, where a different overall factor is used, leading to a sign difference.

tidal_heating(pnsystem)

Compute the rate of energy and angular-momentum absorption into each black hole in a binary.

The returned quantity is a tuple (Ṡ₁, Ṁ₁, Ṡ₂, Ṁ₂), representing the corresponding rates of change of spin (magnitude) and mass on black hole 1 and black hole 2. These apply to general — possibly precessing — non-eccentric binaries. This collection of terms comes from Alvi (2001). It probably wouldn't be too hard to extend Alvi's analysis to eccentric systems.

Note that the validity of the result depends not only on the PN parameter $v$, but also on the angles of the spins relative to the separation vector $n̂$: the smaller the angle, the lower the $v$ at which the approximations should be expected to break down.

See also

• Tagoshi and Sasaki (1994)
• Poisson and Sasaki (1995)
• Tagoshi et al. (1997)
• Porto (2007)
• Chatziioannou et al. (2012)

See the documentation section on "Horizons" for details about the computation of horizon quantities used in this function.

r(pnsystem, [r₀′])
separation(pnsystem, [r₀′])

Compute the separation between the two black holes. This is essentially the multiplicative inverse of γₚₙ, with some factors of G and M thrown in.

r′(pnsystem, [r₀′])
separation_deriv(pnsystem, [r₀′])

Compute the derivative of the separation between the two black holes with respect to v.

Note that we ignore a derivative of M that appears in the r₀′ term, as explained in γₚₙ′.

r⁻¹(r, pnsystem, [r₀′])
separation_inverse(r, pnsystem, [r₀′])

Return v such that r = r(v) when pnsystem is evaluated at v.

Note that the value of v in the input pnsystem is ignored; you may use any value. It may also be convenient to know that you can set the value of v in pnsystem to the returned value using PostNewtonian.vindex as in

pnsystem.state[PostNewtonian.vindex] = r⁻¹(r, pnsystem)

See also γₚₙ⁻¹.

γ̇ₚₙ(pnsystem)
inverse_separation_dot(pnsystem)

Compute the derivative of the inverse separation between the two black holes with respect to time.

γₚₙ(pnsystem, [r₀′])
inverse_separation(pnsystem, [r₀′])

Compute the post-Newtonian parameter

\[\gamma_{\mathrm{PN}} \equiv \frac{G\, M}{r\, c^2},\]

where $r$ is the magnitude of the orbital separation. This quantity has PN order 1, and is given by Eq. (4.3) of Bohé et al. (2013), with spin-squared terms from Eq. (3.32) of Bohé et al. (2015) and spin-cubed terms from Marsat (2014).

Note that there is a 3PN gauge term of $-22ν\ln(r/r₀')/3$. While this value should cancel out of any physical quantity, it is optionally included here for completeness. Computing it requires a few Newton steps to get the value of $γ$ because the $\ln(r)$ term depends on $\gamma$.

Specifically, we use the helper function γₚₙ₀ to write γₚₙ = γₚₙ₀ + (v/c)^8 * (22ln(γₚₙ)/3)ν; given the value of γₚₙ₀, the purpose of this function is to determine γₚₙ.

The default value of r₀′ is 0, in which case that entire term is ignored.

γₚₙ′(pnsystem)
inverse_separation_deriv(pnsystem)

Compute the derivative of γₚₙ with respect to v.

Note that we ignore the time-dependence of the r₀′ term in this function; that constant is obviously independent of v, though it is multiplied by M, which is not independent of v. This dependence, however, should be at a much higher PN order than is currently available in any case, so we ignore it for simplicity.

This computation uses γₚₙ₀ along with the following derivation:

\[\begin{align*} γₚₙ &= γₚₙ₀ + (v/c)^8 (22 \ln(γₚₙ/γ₀′)/3)ν \\ γₚₙ' &= γₚₙ₀' + 8(v/c)^7 (22 \ln(γₚₙ/γ₀′)/3)ν + (v/c)^8 (22 γₚₙ'/3γₚₙ)ν \\ γₚₙ' &= \frac{γₚₙ₀' + 8(v/c)^7 (22 \ln(γₚₙ/γ₀')/3)ν} {1 - (v/c)^8 (22/3γₚₙ)ν} \end{align*}\]

γₚₙ⁻¹(γ, pnsystem, [r₀′])
inverse_separation_inverse(γ, pnsystem, [r₀′])

Return v such that γₚₙ(pnsystem, r₀′) = γ when pnsystem is evaluated at v.

Note that the value of v in the input pnsystem is ignored; you may use any value. It may also be convenient to know that you can set the value of v in pnsystem to the returned value using PostNewtonian.vindex as in

pnsystem.state[PostNewtonian.vindex] = γₚₙ⁻¹(γ, pnsystem)

See also r⁻¹.

γₚₙ₀(pnsystem, [r₀′])

This is a helper function for γₚₙ; this computes the portion of γₚₙ that does not depend on the logarithm of the separation. Since γₚₙ is defined in terms of the separation, that term makes the standard definition into an implicit function.

Specifically, γₚₙ = γₚₙ₀ - (v/c)^8 * (22ln(r*c^2/(G*M)) / 3)ν, where the argument of the logarithm is precisely 1/γₚₙ — by definition. Note that there is a term just like the second term in the expression that involves r₀′, which is a gauge parameter. By default, that term is simply ignored in this function, but if the optional argument is provided, it is included. This combination of the logarithms involving r and r₀′ should drop out of the result for any physical quantity.

γₚₙ₀′(pnsystem)

Helper function to compute γₚₙ′. This just computes the derivative of γₚₙ₀ with respect to v; γₚₙ′ takes care of the extra complications arising from the Newton iterations in γₚₙ′.

ṙ(pnsystem, [r₀′])
separation_dot(pnsystem, [r₀′])

Compute the derivative of the separation between the two black holes with respect to time.

This module contains a few expressions from Kidder (1995).

This is mostly here for testing, because these expressions are not directly used in this package: they are somewhat outdated and describe quantities that are not actually used in this formulation. However, they were used in the SpEC code as an initial guess for eccentricity reduction, so we want to make sure that results from this package are consistent with those from SpEC.

r(pnsystem)

Eq. (4.13).

ṙ(pnsystem)

Eq. (4.12), computed as ṙ = (dE/dt) / (dE/dr), re-expanded and truncated.

aₗ(pnsystem)

Eq. (4.4) of Bohé et al. (2013). This term contributes to 𝛡.

Ω⃗ᵪ(Mⱼ, Mₖ, χ⃗ⱼ, χ⃗ₖ, R)

Compute the angular velocity of precession of spin vector χ⃗ⱼ.

In the approximation that the spin magnitude is constant, the time derivative of $χ⃗ⱼ$ is due to its rotation alone, and is given by $Ω⃗ᵪ × χ⃗ⱼ$.

Note that this function is called by Ω⃗ᵪ₁ and Ω⃗ᵪ₂ with the appropriate parameters; you probably want to use those instead of this one.

The spin-spin term is given by Eq. (2.4) of Kidder (1995); the spin-orbit terms by Eq. (4.5) of Bohé et al. (2013); and the quadrupole-monopole term by Eq. (2.7) Racine (2008).

Ω⃗ᵪ₁(pnsystem)

Compute the angular velocity of precession of χ⃗₁

In the approximation that the spin magnitude is constant, the time derivative of $χ⃗₁$ is due to its rotation alone, and is given by $Ω⃗ᵪ₁ × χ⃗₁$.

Note that this function simply calls Ω⃗ᵪ with the appropriate parameters.

Ω⃗ᵪ₂(pnsystem)

Compute the angular velocity of precession of χ⃗₂

In the approximation that the spin magnitude is constant, the time derivative of $χ⃗₂$ is due to its rotation alone, and is given by $Ω⃗ᵪ₂ × χ⃗₂$.

Note that this function simply calls Ω⃗ᵪ with the appropriate parameters.

Ω⃗ₚ(pnsystem)

Compute the angular velocity of orbital precession.

This is the angular velocity of the orbital angular velocity direction unit vector $ℓ̂$; the time derivative of that unit vector is $Ω⃗ₚ × ℓ̂$.

At the moment, this is computed solely by expressions from Bohé et al. (2013). See 𝛡 for details.

See also R.

𝛡(pnsystem)

Compute the angular velocity of orbital precession according to Bohé et al.

As Bohé et al. (2013) explain above their Eq. (4.1), the orbital precession is given by the time derivative of the orbital axis: 𝓵̇ = 𝛡 × 𝓵, where the angular velocity is along the separation vector 𝓷, so that 𝛡 = ϖ 𝓷. And in turn, they define aₗ ≔ r ω ϖ, where r is the separation and ω is the orbital angular frequency. Then, they define the PN parameter γₚₙ≔M/r and we have Mω = v³ so that ϖ = γₚₙ aₗ / v³. The parameters γₚₙ and aₗ are given by Eqs. (4.3) and (4.4), and given here by the functions γₚₙ and aₗ.

The spin-squared terms (by which we mean both spin-spin and spin-orbit squared terms) in the energy are known to 3pN order, and given in Eq. (3.32) of Bohé et al. (2015).

See also R.

h!(h, pnsystem; ℓₘᵢₙ=0, ℓₘₐₓ=typemax(Int))
mode_weights!(h, pnsystem; ℓₘᵢₙ=0, ℓₘₐₓ=typemax(Int))

Compute mode weights of gravitational waves emitted by pn system, modifying h in place.

These modes are computed in the "co-orbital" frame, in which the larger object lies on the positive $x$ axis, the smaller lies on the negative $x$ axis, and the instantaneous angular velocity is in the positive $z$ direction.

The modes are stored in h in order of increasing $ℓ$ and increasing $m$, with $m$ iterating fastest, all the way up to the highest available mode, $(8,8)$.

Because gravitational waves have spin weight -2, the $(ℓ,m)=(0,0)$, $(1,-1)$, $(1,0)$, and $(1,1)$ modes are always 0. By default, we assume that these modes are nonetheless included in h. If that is not the case, set ℓₘᵢₙ to the smallest $ℓ$ value that should be present in the output data — ℓₘᵢₙ=2 being the most reasonable alternative.

These results come most directly from Eqs. (A5) of Boyle et al. (2014), with the exception of errors in the 2PN spin-spin terms, in which cases we must multiply by $\nu/2$ and make the substitutions $S_1 \mapsto S_0^+$ and $S_2 \mapsto S_0^-$. In turn, those expressions are synthesized from the following: Non-spinning terms are taken from Blanchet (2014), except for the highest-pN term in the (2,±2) mode, which are taken from Blanchet et al. (2023), and the $m=0$ modes, which are taken from Favata (2008). The 1PN spin-orbit term is from Eq. (3.22d) of Kidder (1995). The 1.5PN spin-orbit term is from Eq. (3.22f) of Kidder (1995) and Eq. (F15b) of Will and Wiseman (1996). The 2PN spin-orbit term is from Eq. (4.13) of Buonanno, Faye, Hinderer (2013), while the 2PN spin-spin term is from Eq. (4.15) of that reference.