Add gravitational radiation emission to binary evolution.

[Adam-Boesky]

Add the option to emit gravitational wave emission during binary evolution. To approximate the effects of gravitational radiation, we use Peters 1964 Eq. 5.6 and 5.7 to update the semimajor axis and eccentricity, respectively. We add dynamic timestepping where we set dt = -1e-2 * a / (da/dt) if it is less than the existing dt.

Testing:
For a test BBH binary, we obtain a delay time of 1664.13 Myr compared to the post-processing calculation of 1660.34 Myr for the same binary. This is approximately 0.23 percent error.

Time to evolved 10000 binaries:

• With GWs: 215.692 seconds
• Current COMPAS: 222.634 seconds

[jeffriley]

Thanks Adam! Just a few comments for now:

1. c++ is a strongly typed language. That means when you divide integers, you will get an integer result.

For example, at line 2234 in void BaseBinaryStar::CalculateGravitationalRadiation(), you have:

m_DaDtGW = (numeratorA / denominatorA) * (1 + (73 / 24) * eccentricitySquared + (37 / 96) * eccentricitySquared * eccentricitySquared) * MYR_TO_YEAR; // units of AU Myr^-1

The divisions of the integer constants there might not produce the results you're expecting. Consider the following toy c++ program:

#include <iostream>

int main() {

    double v1 = 73 / 24;
    double v2 = 73.0 / 24.0;
    double v3 = 1.7 * 73 / 24;
    double v4 = 1.7 * (73 / 24);
    double v5 = 1.7 * (73.0 / 24.0);

    std::cout << "73 / 24 = " << v1 << "\n";
    std::cout << "73.0 / 24.0 = " << v2 << "\n";
    std::cout << "1.7 * 73 / 24 = " << v3 << "\n";
    std::cout << "1.7 * (73 / 24) = " << v4 << "\n";
    std::cout << "1.7 * (73.0 / 24.0) = " << v5 << "\n";

    return 0;
}

When that program is compiled, linked, and run, the result is:

73 / 24 = 3
73.0 / 24.0 = 3.04167
1.7 * 73 / 24 = 5.17083
1.7 * (73 / 24) = 5.1
1.7 * (73.0 / 24.0) = 5.17083

It is best to always express floating point constants with the decimal place (and preferably the following 0 - so, 23.0 rather than 23.) so that you can be sure you're actually getting the result you want. Even in expressions like that on lines e.g. 2228, 2232, and 2233. It's just better practice to always put the decimals in - for future readers of the code, and in case someone in the future modifies the expression and gets into trouble because the decimals aren't there and we end up with integer arithmetic that wasn't intended.

2. We quantise the timestep in COMPAS to an integer multiple of the TIMESTEP_QUANTUM (1.0E-12 Myr). You'll need to quantise the result from line 2644 in BaseBinaryStar::Evolve() (see example a few lines above).

3. You need to create a new default yaml file - use './compas --create-yaml-file yourfilename.yaml', and then move 'yourfilename.yaml' to 'compas_python_utils/preprocessing/compasConfigdefault.yaml'
   (see documentation at https://compas.readthedocs.io/en/latest/pages/User%20guide/Pre-processing/pre-processing-yaml.html)

4. Not a big deal, but the parentheses around 'aNew' are not required at line 2254 in BaseBinaryStar::EmitGravitationalWave()

5. You need to add the new option to the documentation (at https://compas.readthedocs.io/en/latest/pages/User%20guide/Program%20options/program-options-list-defaults.html) - note the category listing at the bottom of the page - you'll need to add the new option to the alphabetical list, and the category section at the bottom. Documentation files are reStructuredText (RST) files in the 'online-docs' folder (in this case online-docs/pages/User Guide/Program Options/)

6. You need to document (briefly) what you've added/changed in the code in changelog.h - it's fairly straight-forward. You will need to change the COMPAS version number (at the end of the changelog file). Our convention is M.mm.dd, where M is a major change, mm is a minor change, and dd is a defect repair. In this case, since you are making a minor change to functionality that is not a defect repair, you should increment the middle numbers (mm).

7. You could (I think you should) also add an entry to the "what's-new" page in online-docs/pages just to announce the new option - again, it's fairly straight-forward (just use previous entries as a template).

[ilyamandel]

P.S. @Adam-Boesky -- note that there is a conflict because another PR got merged in the mean time; please merge the latest dev into your branch and resolve any conflicts before resubmitting.

[ilyamandel]

@Adam-Boesky -- thanks for the changes. @jeffriley -- can you check whether these match your requests?

I am surprised Jeff didn't object to the fact that m_DaDtGW and m_DeDtGW were effectively used as mock globals to avoid passing information between two functions; they aren't used anywhere else and don't seem to be loggable -- perhaps they should be? -- but, as it is, I don't see the why they are needed at class-level variables... However, I defer to Jeff.

My one request is about this bit of code:

m_SemiMajorAxis = aNew < 0.0 ? 1E-20 : (aNew); // if <0, set to arbitrarily small number

// Update the eccentricity
m_Eccentricity += m_DeDtGW * p_Dt;

First, what is the reason for setting m_SemiMajorAxis to an arbitrary small number but not exactly zero? Second, there's an asymmetry here: what if eccentricity goes negative? Finally, should we use Jeff's nice comparison function rather than an exact comparison?

[jeffriley]

My first review was just a high-level, first-pass - just to make sure we got the changelog and documentation done (and a couple of things that popped out at me :-)). I'll do a deeper review now that @ilyamandel has had a look.

[jeffriley]

@Adam-Boesky, just a few quick comments in addition to my earlier review and @ilyamandel 's review.

Since other PRs have been merged:
whats-new.rst will need to be updated (version)
changelog.h will need to be updated

I agree with @ilyamandel re m_DaDtGW and m_DeDtGW. BaseBinaryStar::CalculateGravitationalRadiation() should return a tuple containing two doubles rather than use class member variables m_DaDtGW and m_DeDtGW

In BaseBinaryStar::CalculateGravitationalRadiation() I think I'd break the calculations down a bit further and use local variables for things like:

G_AU_Msol_yr_3 = G_AU_Msol_yr * G_AU_Msol_yr * G_AU_Msol_yr;
C_AU_Yr_5 = C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr;

just to make the long calculations a bit more readable and easier for a human to parse.

Similarly for multiple multiplications of m_SemiMajorAxis and oneMinusESq, and the benefit there is that you'd have a common factor - e.g. for m_SemiMajorAxis

m_SemiMajorAxis_3 = m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis;

elsewhere you have m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis, so you could use m_SemiMajorAxis_3 * m_SemiMajorAxis and only do the ^3 once. And similarly for oneMinusESq.

I think the descripion of BaseBinaryStar::EmitGravitationalWave():

/*

• Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation()
• void EmitGravitationalRadiation(const double p_Dt)
  ...

should indicate that the function modifies the sem-major axis, eccentricity, and previous eccentricity values (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev), something like:

/*

• Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation()
• This function updates the sem-major axis, eccentricity, and previous eccentricity values
• (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev) as a result of emitting the GW
• void EmitGravitationalRadiation(const double p_Dt)

and also in BaseBinaryStar::EvaluateBinary():

// emit gravitational wave if required
// this may update m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev

if (OPTIONS->EmitGravitationalRadiation()) {
EmitGravitationalWave(p_Dt);
}

just so people reading the code know that emitting the GW will change the binary (they should, but stating it explicitly will bring it to their attention).

I'll take another look later.

[ilyamandel]

@Adam-Boesky, just pinging you to make sure you've seen the change requests.

[ilyamandel]

BaseBinaryStar.cpp:
Use Peters 1964's approximations of the effects of emitting a GW -> ...emitting GWs

 * void CalculateGravitationalRadiation()
 * 
 * @return                  The change in semi major axis and eccentricity per time due to GW emission

The first line should not say void for the return type; for the second, include units (e.g., "A tuple containing the rate of change in semimajor axis (AU/Myr) and eccentricity (1/Myr) due to GW emission"

 * This function updates the sem-major axis, eccentricity, and previous eccentricity values
 * (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev) as a result of emitting the GW.

sem-major is a typo
Do we need to update m_EccentricityPrev here? I think not: other things can change eccentricity in a given timestep, and don't want to update more than once per timestep.
"emitting GWs"

     m_SemiMajorAxis = aNew < 0.0 ? 1E-20 : (aNew);  // if <0, set to arbitrarily small number

    // Update the eccentricity
    m_Eccentricity += DeDtGW * p_Dt;

See my previous comments.  Should these be caught as mergers?

    
        // Save values as previous timestep
    m_SemiMajorAxisPrev = m_SemiMajorAxis;
    m_EccentricityPrev = m_Eccentricity;

Apparently, you also update m_SemiMajorAxisPrev; see comment above -- don't think this is the right place to do this.

            double DaDtGW;
        double DeDtGW;

See previous comments -- can we avoid these pseudo-globals to pass parameters around?

Line 2895: it seems like this will only correct the timestep on the next step; should we be aborting the current timestep and taking a smaller one if the chosen timestep was too large, or updating dt earlier?

changelog.h:
- If --emit-gravitational-radiation, timestep dynamically as function of da/dt.
Missing verb

Options.cpp:

m_EmitGravitationalRadiation = true;
We should set this to false here, until we have tested things thoroughly and agree to make this the default behavior. (This also requires updating the default in the documentation.)

[jeffriley]

I put the following in the proj_gw_radiation slack channel - repeating it here:

What is it that we want to record about mergers here? The fact? Or details about the merger? We already have a merger flag in the system parameters file - in fact we have two: merger at birth, and merger (after birth). Is what you want to record here just the fact of the merger (so would the merger flag (after birth) flag in the system parameters file suffice, or do you want attributes of the constituent stars and the binary to be recorded at the time of the merger (so we would need a record in (say) the CE file)?

Currently we print a record to the BSE_Common_Envelopes file in BaseBinaryStar::ResolveCommonEnvelopeEvent() when we resolve a common envelope event, and as the code is now we only have one record type for the common envelopes file (the default record type). This PR adds a call to PrintCommonEnvelope() in BaseBinaryStar::Evolve() - meant to record the merger. That call will add a record to the common envelopes file, but users will just interpret it as another CE event - unless we use a different record type for that entry. But see the paragraph above - do we need to print an entry to the common envelopes file, or will flagging the fact that the merger occurred in the system parameters file suffice? If we just want to flag the merger, but we want to differentiate from an ordinary stellar merger, the we could change the Merger files in the system parameters file from a boolean to an integer so that we can differentiate, or we could just add another column for DCO mergers (that's what these are, right?).

So if we just want to record the fact of the merger, remove the added call to PrintCommonEnvelope() in BaseBinaryStar::Evolve() and figure out how we want to flag the merger in the system parameters file. If we really do want to record attributes of the stars and binary when the merger occurs in the common envelopes file, then I think we should add a new record type for the common envelopes file.

[jeffriley]

Line 2895: it seems like this will only correct the timestep on the next step; should we be aborting the current timestep and
taking a smaller one if the chosen timestep was too large, or updating dt earlier?

I vote for updating dt earlier - especially if aborting the timestep means reverting state...

[ilyamandel]

@Adam-Boesky :

To be more specific on setting dt:
I think what we should be doing is updating the guess for dt which is currently made around line 2831 of BaseBinaryStar.cpp (see also the first timestep value around line 2689 of the same file):
dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // calculate new timestep

We should be taking the minimum of the stellar timestep guesses and the guess from the binary's orbital evolution (0.01 of the timescale for the orbit to change). I'd suggest moving this line to a separate function, like CalculateTimestepBinary(), since there will be more than one reason for the binary to evolve (GW radiation, tides, ...) that can all be added to this function.

[jeffriley]

wrt logging, since we're going to log to the CE file, we should add a new record type to the CE file record types - call it DCO_MERGER perhaps? - and add that parameter to your newly-added call to PrintCommonEnvelope(). When adding that, change the default record type name from DEFAULT to COMMON_ENVELOPE (unless there's a better name), and also change DEFAULT to COMMON_ENVELOPE in the declaration of PrintCommonEnvelope() in BaseBinarayStar.h. The record types documentation will also need to be updated to reflect those changes (you'll need to merge with at least v02.48.01 to get the record types documentation).

A couple of other things for now:

In BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW), parameters DaDtGW and DeDtGW should have p_ prefix (so readers know they are parameters and not local variables).

If this staement is going to stay in its current form (see @ilyamandel's comment re arbitrarily small number and merger):

m_SemiMajorAxis = aNew < 0.0 ? 1E-20 : (aNew); // if <0, set to arbitrarily small number

use utils::Compare() for the comparison, and remove the parentheses around aNew (not required).

[ilyamandel]

There's some discussion of merger logging in #660 . But I think dealing with mergers properly may be too much to ask of @Adam-Boesky here. I am actually happy if things are labeled as CE events on mergers without any special treatment for the purposes of this PR, as long as the code doesn't crash, etc. We will return to carefully treating mergers later -- I'm planning to implement mergers in the next few weeks regardless (#209 ).

[jeffriley]

Ok, but my suggestion wasn't really intended to treat the merger correctly, just to differentiate between the CE events we log now, and this new CE event. Do we need to be able to differentiate between them, or do we just want to add an entry to the CE file that looks like any other entry? If the latter, do we need to set the CE details (e.g. m_CEDetails - StellarCEDetailsT) correctly before we log the record?

[ilyamandel]

I think we can log as usual, without a special label. This is a merger. Whether stars are brought to merger by stellar evolution, mass transfer, tides, GWs, or something else is not the thing we need to record in a special way, because it's usually a combination of all of these that are at play.

[jeffriley]

@Adam-Boesky I think you need to merge with the latest dev - I believe the latest version is v02.49.00 which is a few revisions ahead of what your changelog.h is showing.

[ilyamandel]

1. I don't think it's feasible to calculate m_DaDtGW (which is needed to compute the timestep) and calculate the change in semi-major axis and eccentricity (which requires a timestep) in the same function. I believe you need to split this up, compute m_DaDtGW and m_DeDtGW (which only require instantaneous binary properties and don't depend on dt), then determine the correct dt (which may or may not be set by the GW timescale), and only then recompute the new binary properties a' = a + m_DaDtGW * dt, etc.

2. if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || p_Dt < NUCLEAR_MINIMUM_TIMESTEP) {
   p_Dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum
   }

We may want to continue evolving merger products (which would be characterised by one component becoming the merger product and the other being a massless remnant). We don't want to take tiny timesteps in that case; on the contrary, the timestep should be driven exclusively by the dt passed back by the star which is not a massless remnant.

Minor:
The version should be updated in whats-new (02.47.00 May 17, 2024).

[Adam-Boesky]

Responding to @ilyamandel's most recent points:

1. (mentioned in Slack but copied here) The reason that I chose to calculate the m_DaDt/m_DeDt and updated the values in the same step is because I thought it was a sufficient approximation to use the previous dt and it helps avoid passing around variables. I tested this and it had a negligible effect on the delay time of a binary, and it is similar to how we currently handle updating timesteps. I could change it back to a two-function process, but that will require passing around parameters or adding class variables, and the current version is (at least) a good approximation based on preliminary testing.
2. I took this line directly from the current version of BaseBinarStar.cpp. When you say we want to continue evolving merger products, is this specifically for when EmitGravitationalRadiation is True? If so, I could add some sort of condition such that we consider the minimum nuclear timestep only when EmitGravitationalRadiation==False...

[ilyamandel]

Hi @Adam-Boesky,
Sorry, I somehow missed your response from 3 weeks ago and though that you were too busy to respond for now (Floor told me you were working hard for Mr. Musk...). I apologise!
On item 2, I take your point that this is an existing issue, so all good, to be fixed later.
On item 1, I still think it's best to separate out timestep choices and actually taking the timestep. Things may have happened since the last timestep (say, a common envelope or a supernova kick) that would significantly change the binary. Furthermore, there may be other evolutionary pressures (say, tides) that want to dictate their own timesteps, and it's best to have a single place where we compare all of the relevant timestep choices and pick the smallest. So I do think it's useful to split this up. Yes, passing around variables is fine in this case.
If you have some spare cycles, would be nice to close this out so that code doesn't diverge further and create more headaches for you -- but no pressure if you are busy, of course!

[ilyamandel]

Hi @Adam-Boesky ,

This looks very good!

I've made a bunch of small changes directly as commits to your branch -- please check that I haven't broken anything in the process.

Still not sure about the following line: why do we need to keep a artificially positive here, and why don't we do a similar check for e to be between 0 and 1? But I can live with this for now.

m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number

Please merge in dev to make sure changelog is up-to-date, and I'm happy to approve!

[ilyamandel]

@jeffriley -- can you also review this when you get a chance (Adam mentioned he'll be offline starting with next week, so I hope we can finish this quickly)?

[Adam-Boesky]

Two responses to @ilyamandel:

1. The following lines are necessary because of how the prev values are used downstream:

    // Save values as previous timestep	
    m_SemiMajorAxisPrev = m_SemiMajorAxis;	
    m_EccentricityPrev = m_Eccentricity;

2. m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number is necessary because HasStarsTouching() includes a check for the semi-major axis being greater than 0. Therefore without this line, if GWs cause the semi-major axis to be less than zero, binaries may continue evolving instead of stopping.

Thanks for the other edits. Tested and everything looks good to go!

[ilyamandel]

OK, I am happy with this for now, so will approve momentarily, but highlighting three issues for future work:

1. I don't think the following is what we want to do going forward (BaseBinaryStar.cpp):

if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || p_Dt < NUCLEAR_MINIMUM_TIMESTEP) {
p_Dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum
}

We may want to continue evolving merger products (which would be characterised by one component becoming the merger product and the other being a massless remnant). We don't want to take tiny timesteps in that case; on the contrary, the timestep should be driven exclusively by the dt passed back by the star which is not a massless remnant.

2. I take your point about setting semi-major axis to a small but positive number to avoid the binary appearing unbound. In retrospect, I think it was a suboptimal implementation decision to equate a negative semi-major axis with the binary being unbound, and we should fix that in the future.

3. I am worried that we may be resetting m_SemiMajorAxisPrev and m_EccentricityPrev multiple times per time step (e.g., what if we have mass transfer and GW emission happening simultaneously?), so these values could end up being overwritten in a way that may not match the developer's or user's expectations. The longer-term solution is to apply all changes to the orbit (from MT, tides, GWs) once per time step.

[ilyamandel]

OK, approving -- in @jeffriley's court.

[Adam-Boesky]

@jeffriley Thanks for the review! Just squashed your commit, please give it a stamp when you get the chance!

[jeffriley]

@Adam-Boesky @ilyamandel

I have opened PR Adam-Boesky#1

I have not merged it - can you take a look at the comments there and the code changes and merge if you agree with my changes.

[jeffriley]

Looks good - thanks @ilyamandel