NaN Values for 'Mass_Core@CO(SN)' with Latest COMPAS (version 03.07.00)

[lee0517-hub]

Describe the bug
When using the latest version of COMPAS (version 03.07.00), I encountered an issue where the data for 'Mass_Core@CO(SN)' is NaN for some binaries.

Label the issue
urgency_high

To Reproduce
Steps to reproduce the behavior:
./COMPAS --critical-mass-ratio-prescription NONE --kick-magnitude-distribution ZERO --random-seed 1724542090 --common-envelope-alpha 1.6265542767303343 --wolf-rayet-multiplier 0.8084912681292984 --number-of-systems 1 --initial-mass-function KROUPA --mass-ratio-distribution DUQUENNOYMAYOR1991 --semi-major-axis-distribution SANA2012 --minimum-secondary-mass 5 --metallicity-distribution LOGUNIFORM --tides-prescription PERFECT --remnant-mass-prescription FRYER2022 --common-envelope-allow-radiative-envelope-survive False --common-envelope-allow-immediate-RLOF-post-CE-survive False

Expected behavior
A clear and concise description of what you expected to happen.

Screenshots

Versioning (please complete the following information):

• OS: [e.g. Ubuntu 22.04]
• COMPAS [03.07.00]

Additional context
Add any other context about the problem here.

[jeffriley]

The immediate cause of this problem is that in HeGB::CalculateCoreMassOnPhase_Static(), p_time is greater than tx, but it is also greater than tinf2, so we end up taking the sqrt of a negative number, hence the NaN.

The calculation is given by Hurley et al., 2000 eq 39, modified per the discussion following eq 84. I believe our implementation of the equation is correct, so we need to determine why p_Time is > tinf2. Hurley doesn't mention caveats, so I assume he expected tinf2 would always be greater than p_Time (t in eq 39). Preventing the NaN naively is easy, but I'm not sure that would be physical... Needs more investigation.

[lee0517-hub]

It seems that this bug has only recently appeared. I am not certain from which version it started, but I can confirm that the issue does not exist in version v03.00.05. I hope this information helps and saves you some time in the debugging process.

[jeffriley]

Problem was introduced in v03.05.02.

@ilyamandel

Prior to v03.05.02, running the command @lee0517-hub gives above produces a single entry in the BSE_Supernovae file (no NaNs).

Running the same command with v03.05.02 results in a second SN event, with both Mass_Core@CO(SN) and Mass_CO_Core@CO(SN) set to Nan.

I ran the same command with v02.50.01, and there was only a single entry in the BSE_Supernovae file, so I don't believe this problem was introduced by reinstating IsCCSN() etc. to pre-v03.00.00 states.

I removed the call to EvolveOnPhase(0.0) from the Initialise() function of all evolved stellar types in v03.05.02 and re-ran the command shown above, and there was only a single entry in the BSE_Supernovae file, so it would seem the addition of the call to EvolveOnPhase(0.0) in Initialise() for all evolved stellar types has introduced the problem (or maybe just exposed it, but I think we need to determine if the call toEvolveOnPhase(0.0) in Initialise() is producing valid results - possibly making attribute values a step ahead of where they should be?).

EDIT: changelog.h indicates that you added a call to EvolveBinary(0.0) - it was actually EvolveOnPhase(0.0) - we should probably fix changelog.h.

EDIT: @ilyamandel should we be updating timescales and gbparams after calling EvolveOnPhase(0.0) in Initialise()?

[nrsegovia]

Possibly related: long period (>~800d) systems containing and sdB (type 7) and a low mass MS star (Msdb/Mms ~0.4 to 0.8, Msdb ~ 0.47 Msun) have completely vanished after introducing EvolveOnPhase(0.0). I think the relevant case here is the CHeB stage, though I have looked at the detailed evolution and some properties like Age and Radius make more sense after the change: Age does not take negative values anymore, while Radius matches the CHeB size (unlike before, where it followed the FGB). The latter is the reason why these systems existed in the first place, as the RLOF only happened at the very tip of the FGB when COMPAS already flags the star as CHeB. Such RLOF does not happen in the latest version, which is a problem since these systems should exist (they are observed).

One example:

COMPAS --random-seed 15 --maximum-evolution-time 14000 -n 1 --initial-mass-1 1.3699312401840573 --initial-mass-2 0.8546012889119561 --orbital-period 878.3435151466048 --eccentricity 0.3287122621998207 -z 0.014768 --common-envelope-alpha 0.2 --mass-transfer-accretion-efficiency-prescription FIXED --critical-mass-ratio-prescription GE20 --mass-transfer-fa 0.0 --mass-transfer-angular-momentum-loss-prescription MACLEOD_LINEAR --mass-transfer-jloss-macleod-linear-fraction-degen 0 --mass-transfer-jloss-macleod-linear-fraction-non-degen 0 --radial-change-fraction 0.005 --mass-change-fraction 0.005

Before the change, the primary star in this system would follow the stellar type sequence 3-4-7 (RLOF at the end of 3/beginning of 4), now it follows 3-4-5 (no RLOF).

[ilyamandel]

I tried to reproduce @lee0517-hub 's issue, but could not do so either with that command line or with a longer run (no NaNs) when running 03.09.03. It may have been fixed when changing the treatment of mergers in CEs. If it's still an issue with that latest code, I'd be happy to look into this further.

[ilyamandel]

Hi @nrsegovia ,

It seems that the situation you are describing is that we were artificially producing certain systems (albeit ones that should exist) because of buggy code. I am not sure that's a good reason to revert the change.

But it is interesting to understand what the expected evolutionary sequence should be and why we can't reproduce it.

Are you saying that we do form sdB+MS stars, but the masses of the latter are too high? That might point to our needing to assume more stable mass transfer. Or is it that the masses are right, but the separations too small? That may point to less angular momentum loss on non-conservative mass transfer. In other words, I am asking whether this points us to constraints on some evolutionary parameters for which we are currently assuming incorrect values.

Or did I misunderstand you, and are you saying that you think there is a bug in the code (in the sense that the code isn't doing what it was intended to do, rather than doing what was intended but the intention being physically incorrect)?

[nrsegovia]

Are you saying that we do form sdB+MS stars, but the masses of the latter are too high? That might point to our needing to assume more stable mass transfer. Or is it that the masses are right, but the separations too small? That may point to less angular momentum loss on non-conservative mass transfer. In other words, I am asking whether this points us to constraints on some evolutionary parameters for which we are currently assuming incorrect values.

We do expect sdB + massive MS, though they are harder to constrain observationally (the bright MS star would outshine the sdB). Some authors have predicted these systems using other theoretical methods, so I'm not concerned about them showing up in COMPAS. I wouldn't think that angular momentum loss or mass transfer itself is the problem either, the way I see it (after taking a look at the detailed output) is that these systems are rather 'delicate' in terms of time step choices. As the code is now, we get up to R/RL ~ 0.9988 before the star evolves to the CHeB and significantly shrinks its radius. Before, the code would allow an additional time step that would let the star go up to R/RL ~1.0005, triggering a MT episode. This might have been a 'bug' since the star would be labeled as CHeB but still have a Giant-like radius.

I was just thinking that the EvolveOnPhase(0.0) call might have somehow introduced an earlier-than-desired phase change into the CHeB stage, but now I see that it is the other way around and it just updates properties adequately (unless I'm missing something and the phase change should indeed happen a time step later...).

TL;DR: I do not think it is a problem with the physics side of things. It is just a type of system that is really sensitive to time step-related choices. In fact, removing the --radial-change-fraction 0.005 --mass-change-fraction 0.005 part allows RLOF, though the result is a HeWD (type 10) instead.