Dynamics

orbital_evolution(pnsystem; kwargs...)
orbital_evolution(M₁, M₂, χ⃗₁, χ⃗₂, Ωᵢ; kwargs...)

Integrate the orbital dynamics of an inspiraling non-eccentric compact binary.

Required arguments

The first argument to this function may be a single PNSystem that encodes these required arguments (as well as Rᵢ, Λ₁, and Λ₂ among the keyword arguments), or the following may be given explicitly:

• M₁: Initial mass of object 1
• M₂: Initial mass of object 2
• χ⃗₁: Initial dimensionless spin of object 1, S⃗₁/M₁²
• χ⃗₂: Initial dimensionless spin of object 2, S⃗₂/M₂²
• Ωᵢ: Initial orbital angular frequency

(Note that the explicit inputs require Ωᵢ, whereas PNSystems require vᵢ as input.)

These parameters all describe the "initial" conditions. See below for an explanation of the different meanings of "initial" and "first" in this context. Note that the masses change in time as a result of tidal heating — though the changes are quite small throughout most of the inspiral. The spins change direction due to precession, but also change in magnitude due to tidal heating. Therefore, the values passed here are only precisely as given precisely at the moment of the initial data corresponding to the frequency Ωᵢ.

Keyword arguments

Note that several of these keywords are given as Unicode but can also be given as the ASCII string noted. For example, Λ₁ may be input as Lambda1 equivalently; the default values are the same, regardless.

• Λ₁=0 or Lambda1: Tidal-coupling parameter of object 1.
• Λ₂=0 or Lambda2: Tidal-coupling parameter of object 2.
• Ω₁=Ωᵢ or Omega_1: First angular frequency in output data. This may be less than Ωᵢ, in which case we integrate backwards to this point, and combine the backwards and forwards solutions into one seamless output. (See next section.)
• Ωₑ=Ω(v=1,M=M₁+M₂) or Omega_e: Final angular frequency at which to stop ODE integration. Note that integration may stop before the system reaches this frequency, if we detect that PN has broken down irretrievably — for example, if one of the masses is no longer strictly positive, if a spin is super-extremal, or the PN velocity parameter v is decreasing, or is no longer in the range (0,1). Warnings will usually only be issued if v < 0.35, but if quiet=true informational messages will be issued.
• Rᵢ=Rotor(1) or R_i: Initial orientation of binary.
• approximant="TaylorT1": Method of evaluating the right-hand side of the evolution equations.
• PNOrder=typemax(Int): Order to which to retain powers of $v^2$ in PN expansions. The default is to include all available terms in each PN expression.
• check_up_down_instability=true: Warn if the "up-down instability" (see below) is likely to affect this system.
• time_stepper=Vern9(): Choice of solver in OrdinaryDiffEq to integrate ODE.
• abstol=eps(T)^(11//16): Absolute tolerance of ODE solver, where T is the common type to which all the positional arguments are promoted. This is the tolerance on local error estimates, not necessarily the global error. Note that 11//16 is just chosen to suggest that we will have roughly 11 digits of accuracy (locally) for Float64 computations, and a similar accuracy for other float types relative to that type's epsilon.
• reltol=eps(T)^(11//16): Relative tolerance of ODE solver. (As above.)
• termination_criteria_forwards=nothing: Callbacks to use for forwards-in-time evolution. See below for discussion of the default value.
• termination_criteria_backwards=nothing: Callbacks to use for backwards-in-time evolution. See below for discussion of the default value.
• force_dtmin=true: If dt decreases below the integrator's own minimum, and this is false, the integrator will immediately raise an error, before the termination criteria have the chance to exit gracefully. Note that a true value here is critical if the dtmin_terminator callback is to have any effect.
• quiet=true: If set to false, informational messages about successful terminations of the ODE integrations (which occur when the target $v$ is reached in either direction) will be provided. Warnings will still be issued when terminating for other reasons; if you wish to silence them too, you should do something like

  using Logging
  with_logger(SimpleLogger(Logging.Error)) do
      <your code goes here>
  end

• saves_per_orbit=0: If greater than 0, the output will be interpolated so that there are saves_per_orbit time steps in the output for each orbit. Note that this conflicts with the saveat option noted below.

All remaining keyword arguments are passed to the solve function of DiffEqBase. See that function's documentation for details, including useful keyword arguments. The most likely important one is

• saveat: Denotes specific times to save the solution at, during the solving phase — either a time step or a vector of specific times.

In particular, if you want the solution to be output at uniform time steps δt, you want to pass something like saveat=δt; you don't want the solve keyword dt, which is just the initial suggestion for adaptive systems. It is not permitted to pass this option and the saves_per_orbit option.

Also note that callback can be used, and is combined with the callbacks generated by the termination_criteria_* arguments above. That is, you can use the default ones and your own by passing arguments to callback. See the documentation for more details, but note that if you want to make your own callbacks, you will need to add OrdinaryDiffEq to your project — or possibly even DifferentialEquations for some of the fancier built-in callbacks.

ODE system

The evolved variables, in order, are

• M₁: Mass of black hole 1
• M₂: Mass of black hole 2
• χ⃗₁ˣ: $x$ component of dimensionless spin of black hole 1
• χ⃗₁ʸ: $y$ component...
• χ⃗₁ᶻ: $z$ component...
• χ⃗₂ˣ: $x$ component of dimensionless spin of black hole 2
• χ⃗₂ʸ: $y$ component...
• χ⃗₂ᶻ: $z$ component...
• Rʷ: Scalar component of frame rotor
• Rˣ: $x$ component...
• Rʸ: $y$ component...
• Rᶻ: $z$ component...
• v: PN "velocity" parameter related to the total mass $M$ and orbital angular frequency $Ω$ by $v = (M Ω)^{1/3}$
• Φ: Orbital phase given by integrating $Ω$

The masses and spin magnitudes evolve according to tidal_heating. The spin directions evolve according to Ω⃗ᵪ₁ and Ω⃗ᵪ₂. The frame precesses with angular velocity Ω⃗ₚ, while also rotating with angular frequency Ω about the Newtonian orbital angular velocity direction. The frame rotor $R$ is given by integrating the sum of these angular velocities as described in Boyle (2016). And finally, the PN parameter $v$ evolves according to something like

\[\dot{v} = - \frac{\mathcal{F} + \dot{M}_1 + \dot{M}_2} {\mathcal{E}'}\]

where 𝓕 is the flux of gravitational-wave energy out of the system, $\dot{M}_1$ and $\dot{M}_2$ are due to tidal coupling as computed by tidal_heating, and 𝓔′ is the derivative of the binding energy with respect to $v$. For "TaylorT1", the right-hand side of this equation is evaluated as given; for "TaylorT4", the right-hand side is first expanded as a Taylor series in $v$ and then truncated at some desired order; for "TaylorT5", the inverse of the right-hand side is expanded as a Taylor series in $v$, truncated at some desired order, and then inverted to obtain an expression in terms of $v$.

Returned solution

The returned quantity is an ODESolution object, which has various features for extracting and interpolating the data. We'll call this object sol.

The field sol.t is the set of time points at which the solution is given. To access the ith variable at time step j, use sol[i, j].^[1] You can also use colons. For example, sol[:, j] is a vector of all the data at time step j, and sol[i, :] is a vector of the ith variable at all times.

For convenience, you can also access the individual variables with their symbols. For example, sol[:v] returns a vector of the PN velocity parameter at each time step. Note the colon in :v, which is Julia's notation for a Symbol.

Initial frequency vs. first frequency vs. end frequency

Note the distinction between Ωᵢ (with subscript i) and Ω₁ (with subscript 1). The first, Ωᵢ, represents the angular frequency of the initial condition from which the ODE integrator will begin; the second, Ω₁, represents the target angular frequency of the first element of the output data. That is, the ODE integration will run forwards in time from Ωᵢ to the merger, and then — if Ωᵢ>Ω₁ — come back to Ωᵢ and run backwards in time to Ω₁. The output data will stitch these two together to be one continuous (forwards-in-time) data series.

For example, if you are trying to match to a numerical relativity (NR) simulation, you can read the masses and spins off of the NR data when the system is orbiting at angular frequency Ωᵢ. Integrating the post-Newtonian (PN) solution forwards in time from this point will allow you to compare the PN and NR waveforms. However, you may want to know what the waveform was at earlier times than are present in the NR data. For this, you also have to integrate backwards in time. We parameterize the point to which you integrate backwards with Ω₁. In either case, element 1 of the output solution will have frequency Ω₁ — though by default it is equal to Ωᵢ.

Similarly, the optional argument Ωₑ=1 is the frequency of the end element of the solution — that is Julia's notation for the last element. Note that this is automatically reduced if necessary so that the corresponding PN parameter $v$ is no greater than 1, which may be the case whenever the total mass is greater than 1.

Up-down instability

Be aware that the up-down instability (where the more massive black hole has spin aligned with the orbital angular velocity, and the less massive has spin anti-aligned) can cause systems with nearly zero precession at the initial time to evolve into a highly precessing system either at earlier or later times. This is a real physical result, rather than a numerical issue. If you want to simulate a truly non-precessing system, you should explicitly set the in-place components of spin to precisely 0. By default, we check for this condition, and will issue a warning if it is likely to be encountered for systems with low initial precession. The function used to compute the unstable region is up_down_instability.

Time-stepper algorithms

Tsit5() is a good default choice for time stepper when using Float64 with medium-low tolerance. If stiffness seems to be impacting the results, AutoTsit5(Rosenbrock23()) will automatically switch when stiffness occurs. For tighter tolerances, especially when using Double64s, Vern9() or AutoVern9(Rodas5P()) are good choices. For very loose tolerances, as when using Float32s, it might be better to use OwrenZen3().

Termination criteria

The termination criteria are vital to efficiency of the integration and correctness of the solution. The default values for forwards- and backwards-in-time evolution, respectively, are

CallbackSet(
    termination_forwards(v(Ω=Ωₑ, M=M₁+M₂)),
    dtmin_terminator(T),
    decreasing_v_terminator(),
    nonfinite_terminator()
)

and

CallbackSet(
    termination_backwards(v(Ω=Ω₁, M=M₁+M₂)),
    dtmin_terminator(T),
    nonfinite_terminator()
)

where T is the common float type of the input arguments. If any additional termination criteria are needed, they could be added as additional elements of the CallbackSets. See the callback documentation for details.

uniform_in_phase(solution, saves_per_orbit)

Interpolate solution to uniform steps in phase.

By default, the solution returned by orbital_evolution may be sampled very sparsely — too sparsely to satisfy the Nyquist limit of the waveform. If the waveform extends to $\ell_{\mathrm{max}}$, there will be modes varying slightly more rapidly than $\exp\left(\pm i\, \ell_{\mathrm{max}}\, \Phi \right)$, where $\Phi$ is the orbital phase. If the frequency were constant, this would require at least $2\ell_{\mathrm{max}}$ samples per orbit. To incorporate a safety factor, $4\ell_{\mathrm{max}}$ seems to work fairly reliably.

See also the saves_per_orbit and saveat arguments to orbital_evolution, as well as interpolation-in-time capabilities of the result of that function.

estimated_time_to_merger(M, ν, v)
estimated_time_to_merger(pnsystem)

Compute the lowest-order PN approximation for the time to merger, starting from PN velocity parameter v.

This is used internally as a convenient way to estimate how long the inspiral integration should run for; we don't want it to integrate forever if PN has broken down. However, it can be a very poor approximation, especially close to merger, and doubly so if the spins or eccentricity are significant.

fISCO(q, M)
fISCO(pnsystem)

Compute the "BKL" approximation for the ISCO (Innermost Stable Circular Orbit) frequency.

This is taken from Eq. (5) of Hanna et al. (2008). Note that this does not account for the spins of the objects in the binary, so that this returns a very crude estimate of a frequency of interest.

ΩISCO(q,M)
ΩISCO(pnsystem)

2π times fISCO.

up_down_instability(pnsystem)

Compute the range of frequencies over which the system is unstable to increasing precession.

The returned value is a pair of dimensionless frequencies giving the lower and upper frequencies between which we can expect instability. If there is no instability expected, the returned pair is just $(1/M, 1/M)$ — where $M$ is the total mass of the system, and $1/M$ is the upper limit of physically reasonable frequencies.

For compact binaries in which the spins are either aligned or anti-aligned with the orbital angular velocity, we do not expect any precession effects — simply by symmetry. However, if the spin of the higher-mass object is aligned with the orbital angular velocity and the spin of the lower-mass object is anti-aligned, the binary is unstable to precession — meaning that any minuscule misalignment can grow rapidly into significant precession. This was first reported by Gerosa et al. (2015), and the range over which the system is unstable is given by Eq. (2) of that reference. We use the lowest-order approximation to convert binary separation to frequency. The result is also "clamped" between $0$ and $1/M$, because sometimes the PN approximations involved break down and return values outside of those physically plausible limits.

Note that Gerosa et al. use the convention that $q = M_2/M_1$ — which is opposite to the convention used in this package; which we account for internally in this function. They also assume that $M_1 \geq M_2$, which we deal with by automatically swapping the relevant quantities. Neither of these requires any adjustment by users of this function.

function up_down_instability_warn(pnsystem, v₁, vₑ, vₗᵢₘᵢₜ=1//2)

If this system is likely to encounter the up-down instability, log a warning with details.

This function issues the warning if the system is reasonably non-precessing ($\chi_\perp \leq 10^{-2}$) in its current configuration (as given by pnsystem), and the range of frequencies (v₁, vₑ) over which it will be integrated is likely to encounter the up-down instability — except that frequencies above vₗᵢₘᵢₜ will be ignored, as PN is likely to have broken down anyway.

See up_down_instability for details of the calculation of the unstable region.

TaylorT1!(u̇, pnsystem)

Compute the right-hand side for the orbital evolution of a non-eccentric binary in the "TaylorT1" approximant.

This approximant is the simplest, in which the time derivative $\dot{v}$ is given directly by

\[\dot{v} = -\frac{\mathcal{F} + \dot{M}_1 + \dot{M}_2} {\mathcal{E}'},\]

and the PN expression for each term on the right-hand side is evaluated numerically before insertion directly in this expression. Compare TaylorT4! and TaylorT5!.

Here, u̇ is the time-derivative of the state vector, which is stored in the PNSystem object p.

TaylorT4!(u̇, pnsystem)

Compute the right-hand side for the orbital evolution of a non-eccentric binary in the "TaylorT4" approximant.

In this approximant, we compute $\dot{v}$ by expanding the right-hand side of

\[\dot{v} = -\frac{\mathcal{F} + \dot{M}_1 + \dot{M}_2} {\mathcal{E}'}\]

as a series in $v$, truncating again at the specified PN order, and only then is the result evaluated. Compare TaylorT1! and TaylorT5!.

Here, u is the ODE state vector, which should just refer to the state vector stored in the PNSystem object p. The parameter t represents the time, and will surely always be unused in this package, but is part of the DifferentialEquations API.

TaylorT5!(u̇, pnsystem)

Compute the right-hand side for the orbital evolution of a non-eccentric binary in the "TaylorT5" approximant.

In this approximant, we compute $\dot{v}$ by expanding the right-hand side of the inverse of the usual expression

\[\frac{1}{\dot{v}}=\frac{dt}{dv} = -\frac{\mathcal{E}'} {\mathcal{F} + \dot{M}_1 + \dot{M}_2}\]

as a series in $v$, truncating again at the specified PN order, evaluating the result, and then taking the inverse. This approximant was introduced by Ajith (2011). Compare TaylorT1! and TaylorT5!.

Here, u is the ODE state vector, which should just refer to the state vector stored in the PNSystem object p. The parameter t represents the time, and will surely always be unused in this package, but is part of the DifferentialEquations API.

causes_domain_error!(u̇, p)

Ensure that these parameters correspond to a physically valid set of PN parameters.

If the parameters are not valid, this function should modify u̇ to indicate that the current step is invalid. This is done by filling u̇ with NaNs, which will be detected by the ODE solver and cause it to try a different (smaller) step size.

Currently, the only check that is done is to test that these parameters result in a PN parameter v>0. In the future, this function may be expanded to include other checks.