[feat] CO Peters equation integration variable change

[maxbriel]

We use Peters' equations to find if a system merges in the Hubble time or find their final orbital configuration.

Due to the large change required in orbit and time, we need to solve the equations several times to "converge" on our wanted state (for example, contact between Schwarzschild radii for BHs). These multiple steps messes with the other parts of our code, where we set the time and orbital configuration at different moments in DCO evolution.

A change to dimensionless variables only partly solve this, because we're still left with a $1/\alpha^3$ and $1/\alpha^4$ dependencies in the new ODEs, where $\alpha=a/a_0$ (See notes). Thus, as $\alpha \rightarrow 0$ these terms blow up and a very small $\tau$ is required ($\tau=t/t_0$), essentially requiring one below machine precision.

A full change of integration variable is required to get closer to the moment of inspiral. This can be achieved by substituting $s=- \ln(\alpha)$. We want to get close to the moment of merger (contact), which means we never reach $a =0$ but only get as close as possible. As such, we can rewrite our ODEs using this to get close to the final merger and have the solver converge correctly.

Original Peters equations

$$\left< \frac{de}{dt}\right> = -\frac{304}{15}\frac{G^3m_1m_2(m_1+m_2)}{c^5a^4(1-e^2)^{5/2}} e\left(1+\frac{121}{304}e^2 \right)$$

$$ \left<\frac{da}{dt}\right> = -\frac{64}{5} \frac{G^3 m_1 m_2 (m_1 + m_2)}{c^5 a^3 \left(1 - e^2\right)^{7/2}} \left(1 + \frac{73}{24} e^2 + \frac{37}{96} e^4\right) $$

Dimensionless Equations

These are nearly the same as Andrews+19 with a factor 1/4 not embedded in $t_0$

• $\alpha = a/a_0$
• $\tau = t/t_0$
• $t_0 = \frac{5}{64}\frac{c^5a_0^4}{G^3(m_1+m2)m_1m_2} \rightarrow \frac{d\tau}{dt} = \frac{64}{5}\frac{G^3(m_1+m_2)m_1m_2}{c^5a_0^4}$

$$ \frac{d\alpha}{d\tau} = \frac{1}{a_0}\frac{da}{dt}\frac{dt}{d\tau} = - \frac{1}{\alpha^3} \frac{(1+\frac{73}{24}e^2 +\frac{37}{96}e^4)}{(1-e^2)^{7/2}} $$

$$ \frac{de}{d\tau} = \frac{de}{dt} \frac{dt}{d\tau} = -\frac{19}{12}\frac{1}{\alpha^4}\frac{e(1+\frac{121}{304}e^2)}{(1-e^2)^{5/2}} $$

Substituted Equations

• $s=-\ln(\alpha)$
• $\frac{ds}{d\alpha} = -\frac{1}{\alpha}$

$$ \frac{d\tau}{ds} = -\alpha\frac{d\tau}{d\alpha} = e^{-4s} \frac{(1-e^2)^{7/2}}{1+\frac{73}{24}e^2+\frac{37}{96}e^4} $$

$$ \frac{de}{ds} = \frac{de}{d\tau}\frac{d\tau}{ds} = -\frac{19}{12}e(1-e^2)\frac{(1+\frac{121}{304}e^4)}{(1+\frac{73}{24}e^2+\frac{37}{96}e^4} $$

This makes the $de/ds$ independent of $s$ as well.

[astro-abhishek]

We can consider an alternative choice: the analytical fits to Peters given by Mandel 2021: https://ui.adsabs.harvard.edu/abs/2021RNAAS...5..223M/abstract

The accuracy is within 3% (see Fig. 1).

[maxbriel]

We can't use that here sadly, because we have to find the properties (orbit) of living double CO systems as well. For example, LISA populations of double WD or BNS systems.

We also need to have an "interpolated" output to be able to request specific moments in time (in update_after_evo), which also connects to the n_o_steps_history output parameter.

[maxbriel]

Here are the test binaries
evolve_binaries_mb_unitless.txt

And v2.3 from today for comparison.
evolve_binaries_v2.3.txt

There are some differences in the time calculation of systems that reach CO_contact. Most are of the range of ~10 years, though there's one that is a bigger.

[sgossage]

It would be good to have this mirror the other detached_step procedures as much as possible I think.

Here, I'd suggest writing self.evo in such a way that you override the __call__ method of that class for example, so that we just pass and call

    res = solve_ivp(self.evo,
                             ...

here, as in detached_step and its other children classes. Along these lines, I suggest constructing double_CO_evolution to mirror the format of the parent class as well. In the __call__ method, check

    if do_gravitational_radiation:
        # call rhs() here
        ...

    return result

I think this will help maintain modularity and a cohesive format, for example, if others intend to create variations of the detached_step.