TeamCOMPAS / COMPAS

COMPAS rapid binary population synthesis code
http://compas.science
MIT License
64 stars 66 forks source link

BHs are losing mass #884

Closed reinhold-willcox closed 1 year ago

reinhold-willcox commented 1 year ago

Describe the bug BHs appear to be losing mass following the Vink(?) prescription. This was discovered by Andreas Konstantinou. The m_Mdot parameter appears to have related, anomalous behavior.

Label the issue urgency_moderate - This is a moderately urgent issue severity_minor - This is a minor bug with minimal impact

To Reproduce

COMPAS --random-seed 1 --number-of-systems 10000 --metallicity 0.0142 --kick-magnitude-distribution MULLERMANDEL --remnant-mass-prescription  MULLERMANDEL --muller-mandel-kick-multiplier-NS 520 --logfile-definitions COMPAS_Output_Definitions.txt

Together with the following COMPAS_Output_Definitions.txt

BSE_RLOF_Rec += { BINARY_PROPERTY::SYSTEMIC_SPEED, 
    STAR_1_PROPERTY::LUMINOSITY,
    STAR_2_PROPERTY::LUMINOSITY,
    STAR_1_PROPERTY::TEMPERATURE,
    STAR_2_PROPERTY::TEMPERATURE,
    STAR_1_PROPERTY::ANGULAR_MOMENTUM,
    STAR_2_PROPERTY::ANGULAR_MOMENTUM,
    STAR_1_PROPERTY::DYNAMICAL_TIMESCALE,
    STAR_2_PROPERTY::DYNAMICAL_TIMESCALE,
    STAR_1_PROPERTY::THERMAL_TIMESCALE,
    STAR_2_PROPERTY::THERMAL_TIMESCALE,
    STAR_1_PROPERTY::NUCLEAR_TIMESCALE,
    STAR_2_PROPERTY::NUCLEAR_TIMESCALE,
    STAR_1_PROPERTY::MDOT,
    STAR_2_PROPERTY::MDOT,
    STAR_1_PROPERTY::MASS_TRANSFER_DIFF,
    STAR_2_PROPERTY::MASS_TRANSFER_DIFF}

Then mask for accreting BHs.

Expected behavior BHs should not lose mass, and m_Mdot should accurately reflect the mass loss in a given timestep.

Screenshots image

Versioning (please complete the following information):

ilyamandel commented 1 year ago

See Remnants::CalculateMassLossRateHurley(); this is not over-written for NSs or BHs. Hurley+ aren't very clear on whether they intended for this to apply to WDs, but it certainly shouldn't apply to NSs or BHs, so should be overwritten to zero there. Even if Nieuwenhuijzen & de Jager winds are used, they would be zero for BHs because they scale with luminosity -- but a related problem is that our BHs inherit the non-zero luminosity of NSs, which should also be fixed. @akonstantinu would you like to try making these fixes? [@SimonStevenson , pinging you too, since you are looking at winds now.]

ilyamandel commented 1 year ago

@akonstantinu, I still don't understand why the mass loss rates are so high given the expression which I think we are using to compute the Nieuwenhuijzen & De Jager winds -- please check. Also, while you are at it, can you fix NS::CalculateLuminosityOnPhase_Static() -- it's missing a power of 2/3 on the mass relative to Eq. (93) of Hurley, Pols, Tout (2000), which it's supposed to implement.

akonstantinu commented 1 year ago

Do I just make the changes on my PC and push them?

ilyamandel commented 1 year ago

Andreas, you will want to create a branch, test it, create a pull request from that into dev, and after that's reviewed, it will be incorporated into the main code. @reinhold-willcox had a very nice description of the procedure... but now I can't find it, perhaps because it's getting to close to midnight. :) Maybe he can help!

reinhold-willcox commented 1 year ago

Hi @akonstantinu,

Here is the (slightly outdated) workflow. It should get you up to speed on using git and setting up a pull request. If you have any questions, you can ping me.

https://compas.readthedocs.io/en/latest/pages/Getting%20started/dev-git-workflow.html

ilyamandel commented 1 year ago

Working on this issue now. Fixed Nieuwenhuijzen & De Jager winds -- Chris Tout confirms they shouldn't apply to any remnants -- and the missing power in NS::CalculateLuminosityOnPhase_Static().

ilyamandel commented 1 year ago

OK, three things:

  1. Our BHs aren't actually losing mass (at least with the fixes above).
  2. For single stars, Mdot was incorrectly reported for BHs because of the line
    void ResolveMassLoss() { } // NO-OP in Remnants.h. It meant we never called ResolveMassLoss for these, so Mdot was not being reset correctly, and just held onto its last value.
  3. For binary stars, the problem was with BaseBinaryStar::CalculateWindsMassLoss(), where we stopped winds for binaries in mass transfer — but never updated Mdot
ilyamandel commented 1 year ago

Resolved in commit #900 , closing.