Updates to WD accretion to address issue #1384

[SimonStevenson]

Updates to WD accretion to address issue #1384

This enables AIC for ONeWD

Test binary:

./COMPAS -n 1 --mode BSE --initial-mass-1 7.3 --initial-mass-2 5.7 --semi-major-axis 9.0 --evolve-double-white-dwarfs TRUE --emit-gravitational-radiation TRUE --detailed-output TRUE --evolve-unbound-systems FALSE --maximum-neutron-star-mass 2.3

@nrsegovia did you want to add anything else to this PR?

[SimonStevenson]

Running a population test with this PR v dev to test the numbers of SN_Type(SN) == 32 (AIC) and 64 (SN1A):

./COMPAS -n 50000 --mode BSE --evolve-double-white-dwarfs TRUE --emit-gravitational-radiation TRUE --evolve-unbound-systems FALSE --maximum-neutron-star-mass 2.3

AIC:
this PR: 68
dev: 0

SN1a:
this PR: 76
dev: 0

[nrsegovia]

@SimonStevenson I have tweaked the documentation a bit and also added the qCrit fix for one of the prescriptions. I did not remember the variable naming style when I shared the code for the new radius function, so I also modified that. I think it should be all good, though perhaps you might want to add a word or two into ONeWD.cpp explaining that we have assumed that they behave the same way as CO WDs when accreting mass.

[SimonStevenson]

Thanks for the additions Nicolas! Does your fix to qCrit resolve #1385 ?

I do have a comment in the documentation for ONeWD::CalculateMassAcceptanceRate saying

* This currently uses the same prescription as for COWDs, but we may consider different 
 * prescriptions in the future.

I think that should be sufficient for now.

[nrsegovia]

Whoops, I saw the same text as in COWD.cpp and assumed nothing had been specified about it. My bad.

And no, the change does not address the GE prescription from #1385, only HURLEY_HJELLMING_WEBBINK (I implemented this and the Remnant-specific values were removed in one of my cleanups).

[ilyamandel]

Thank you, @SimonStevenson and @nrsegovia !

Comments:

WhiteDwarfs.cpp:

double preSecondFactor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass;

Change 1 to 1.0

---

Remnants.h

double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.59; }

Please add reference for the 1.59 value as a comment.
Is this meant to apply to all Remnants (including NS and BH), or just WhiteDwarfs?

---

If ONeWD.cpp uses the same code as COWD.cpp, we should just call COWD::CaclulateMassAcceptanceRate() from ONeWD.h -- no need to repeat code.

[jeffriley]

Just a comment. We have this code in ONeWD.cpp:

DBL_DBL ONeWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) {

    m_AccretionRegime = DetermineAccretionRegime(p_IsHeRich, p_DonorMassRate); 

Why are those parameters in the reverse order in those functions? I know it doesn't affect the functionality, but I think consistency results in fewer mistakes. Just saying...

[jeffriley]

Cheeky request :-) As long as we're making changes to the WD files, can we also do issue #1351 ?

[nrsegovia]

WhiteDwarfs.cpp:

double preSecondFactor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass;

Change 1 to 1.0

Done, also moved renamed MP and moved it to constants.h

Remnants.h

double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.59; }

Please add reference for the 1.59 value as a comment. Is this meant to apply to all Remnants (including NS and BH), or just WhiteDwarfs?

I have moved all these qCrits to constants.h, perhaps that is better? About remnants: I would think that it only applies to WDs, though the original BSE code does not seem to properly cover other remnants: the else statement in the relevant code block is meant for main sequence stars and sets qCrit to 3. I think that when I first implemented this I also took a look at other codes and, at the very least, COSMIC sets the constant for all remnants. We can default to the 0.0 value as done for CalculateCriticalMassRatio, though. In that case, I would move the relevant line to WhiteDwarfs.h (assuming that no other changes are required).

If ONeWD.cpp uses the same code as COWD.cpp, we should just call COWD::CaclulateMassAcceptanceRate() from ONeWD.h -- no need to repeat code.
I think this should be fine (I have not implemented it).

I will try to address Jeff's comments during the week

[ilyamandel]

We can default to the 0.0 value as done for CalculateCriticalMassRatio, though. In that case, I would move the relevant line to WhiteDwarfs.h (assuming that no other changes are required).

Indeed, I don't think it makes sense to model mass transfer from NS or BH in COMPAS, so we should always label such systems as mergers, regardless of the mass ratio.

I'll wait until you are done with all of the changes you wanted to make before signing off on the review.

[SimonStevenson]

Making ONeWD::CalculateMassAccretionRate call COWD::CalculateMassAccretionRate seems trickier than it may seem. I tried makingCOWD::CalculateMassAccretionRate public and static, but COWD::CalculateMassAccretionRate calls several functions in WhiteDwarfs and I get lots of non-static errors:

COWD.cpp:22:5: error: invalid use of member 'm_AccretionRegime' in static member function
    m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich);
    ^~~~~~~~~~~~~~~~~
COWD.cpp:22:25: error: call to non-static member function without an object argument
    m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich);
                        ^~~~~~~~~~~~~~~~~~~~~~~~
COWD.cpp:27:40: error: call to non-static member function without an object argument
    acceptanceRate = p_DonorMassRate * CalculateEtaHe(p_DonorMassRate);
                                       ^~~~~~~~~~~~~~
COWD.cpp:28:40: error: call to non-static member function without an object argument
    if (!p_IsHeRich) acceptanceRate *= CalculateEtaH(p_DonorMassRate);

I didn't figure out how to resolve this yet. Those functions call member variables so can't be made static.

[nrsegovia]

Indeed, I don't think it makes sense to model mass transfer from NS or BH in COMPAS, so we should always label such systems as mergers, regardless of the mass ratio.

I have moved the qCrits to WhiteDwarfs.h, I believe that now it should fallback to 0.0 for other remnants.

Why are those parameters in the reverse order in those functions? I know it doesn't affect the functionality, but I think consistency results in fewer mistakes. Just saying...

I decided to flip DetermineAccretionRegime parameters. I think all related code has been properly updated, but please take a look.

Cheeky request :-) As long as we're making changes to the WD files, can we also do issue #1351 ?

I found constants in some files and moved them. Let me know if you think I missed any.

I didn't figure out how to resolve this yet. Those functions call member variables so can't be made static.

Perhaps we could move DetermineAccretionRegime to WhiteDwarfs and inherit from there? As long as it can be overloaded in the HeWD case

[SimonStevenson]

Perhaps we could move DetermineAccretionRegime to WhiteDwarfs and inherit from there? As long as it can be overloaded in the HeWD case

Thanks, that seems to work!

[SimonStevenson]

I've merged in dev and updated the changelog and what's new pages. I think this should be ready to go now?

[ilyamandel]

I have moved the qCrits to WhiteDwarfs.h, I believe that now it should fallback to 0.0 for other remnants.

I don't think so, because NeutronStars inherit from Remnants which inherit from He stars. You probably want to explicitly include something like
double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.0; }
in NeutronStar.h.

On identical CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) for ONeWD and COWD: this is declared as a virtual function in WhiteDwarfs.h, so there's absolutely nothing to do for ONeWD -- don't even mention it in ONeWD.h or ONeWD.cpp, it will automatically inherit the COWD version.

I don't understand what you've done with DetermineAccretionRegime(). HeWD inherit from WhiteDwarfs and over-writes it, then COWD inherits from HeWD and doesn't over-write it? I don't think that's the intended behaviour, and I think you do want separate DetermineAccretionRegime() code in COWD if it's different from the HeWD version. The HeWD version, meanwhile, can sit either there or in WhiteDwarfs.

Looks good otherwise, so will approve once these three changes are implemented.

[nrsegovia]

I don't think so, because NeutronStars inherit from Remnants which inherit from He stars. You probably want to explicitly include something like double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.0; } in NeutronStar.h.

Right, I missed the inheritance from HeGB.

On identical CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) for ONeWD and COWD: this is declared as a virtual function in WhiteDwarfs.h, so there's absolutely nothing to do for ONeWD -- don't even mention it in ONeWD.h or ONeWD.cpp, it will automatically inherit the COWD version.

I don't understand what you've done with DetermineAccretionRegime(). HeWD inherit from WhiteDwarfs and over-writes it, then COWD inherits from HeWD and doesn't over-write it? I don't think that's the intended behaviour, and I think you do want separate DetermineAccretionRegime() code in COWD if it's different from the HeWD version. The HeWD version, meanwhile, can sit either there or in WhiteDwarfs.

What I had in mind was to define DetermineAccretionRegime() in WhiteDwarfs, which would then be used by both CO and ONe WDs. As for HeWDs, they have a different behavior so the function would be overloaded. Isn't this what has been implemented? I do not see the inheritance from HeWDs, though I'm no C++ expert. I have the same question for the COWD-ONeWD inheritance in the other comment, particularly since Fig. 1 in the first methods paper does not seem to imply such relation.

[SimonStevenson]

Hi Ilya, thanks for taking a look.

I don't understand what you've done

Honestly I'm not sure I know what I did either, but I don't think it was what I intended to do. I have now edited the code to what I originally intended to do (whether it is correct or not remains to be seen).

ONeWD.h or ONeWD.cpp, it will automatically inherit the COWD version.

ONeWD does not inherit from COWD, it inherits from WhiteDwarfs

then COWD inherits from HeWD

COWD does not inherit from HeWD, it inherits from WhiteDwarfs

The inheritance is not:

WhiteDwarfs -> HeWD -> COWD -> ONeWD

it is:

WhiteDwarfs -> HeWD
WhiteDwarfs -> COWD
WhiteDwarfs -> ONeWD

Therefore, the logic is the following:

DetermineAccretionRegime, CalculateMassAcceptanceRate, CalculateEtaH and CalculateEtaHe are defined in WhiteDwarfs.h and implemented in WhiteDwarfs.cpp.

They are used directly (via inheritance) for COWD and ONeWD.

DetermineAccretionRegime and CalculateMassAcceptanceRate are redefined/overloaded in HeWD which uses a different prescription.

Whilst it may seem counterintuitive that the default implementation is for COWD/ONeWD and not HeWD, this is equally valid given our inheritance order, and avoids the need for duplicate code.

Given the significant changes, I have rerun my tests.

The test binary evolves the same as in my plots in #1384

The population statistics for AIC and SN1A are the same as reported above.

[ilyamandel]

Dear @SimonStevenson , @nrsegovia : you are both correct that COWDs and ONeWDs inherit from WhiteDwarfs, not from HeWDs. I was wrong, of course. :( My sincere apologies. Let me check again.

[ilyamandel]

Looks very good!

Just one question. WhiteDwarfs.h now includes the following comment next to the declaration of DetermineAccretionRegime():
"Can also change flags related to SN events"
Does this refer only to m_HeShellDetonation, or some other flags as well? Might be useful to spell this out.

Other than my optional suggestion above, this is an approval -- happy to see this merged in once @jeffriley OKs it.

[ilyamandel]

Oh, one more thing: What happens if a binary with a WD donor does not satisfy the HjellmingWebbink q-crit threshold? It goes into a CE, right? Do we actually have a meaningful CE implementation for WD donors?

[nrsegovia]

Oh, one more thing: What happens if a binary with a WD donor does not satisfy the HjellmingWebbink q-crit threshold? It goes into a CE, right? Do we actually have a meaningful CE implementation for WD donors?

Good point, I would argue that instead of going into a CE we could assume a merger (from Marsh+ 2004 "Mass transfer is dynamically unstable and will lead to a merger, regardless of synchronization torques, if ... q >2/3").

[ilyamandel]

I agree, @nrsegovia , but is that what would happen now in the code?

[nrsegovia]

@ilyamandel I have now added a couple of lines that should deal with this, though I'm not entirely sure whether this is the best location in the code to do so, and I also wonder how this approach interacts with ResolveCommonEnvelopeEvent

[ilyamandel]

Thank you, @nrsegovia ! But I think your point above, which I agree with, is more general than this particular prescription: any CE initiated by a white dwarf donor should lead to a merger. So I propose the following instead, in the BaseBinaryStar::ResolveCommonEnvelopeEvent() function:

change from

if (isDonorMS || (!envelopeFlag1 && !envelopeFlag2)) {                                                              // stellar merger
    m_MassTransferTrackerHistory = MT_TRACKING::MERGER; 
    m_Flags.stellarMerger        = true;
}

to

if (isDonorMS || (!envelopeFlag1 && !envelopeFlag2) || (m_Star1->IsRLOF() && m_Star1->IsOneOf(WHITE_DWARFS)) || (m_Star2->IsRLOF() && m_Star2->IsOneOf(WHITE_DWARFS)) {                                                              // stellar merger
    m_MassTransferTrackerHistory = MT_TRACKING::MERGER; 
    m_Flags.stellarMerger        = true;
}

[nrsegovia]

@ilyamandel right, makes sense. I'm not familiar with how the other qCrit prescriptions deal with WDs, but I agree that it is better to have a "standardized" behavior. I have added the code you suggested and removed the couple of lines I added before.

[ilyamandel]

Thank you, @nrsegovia !

I'll approve, but you need to merge the current version of dev into this (it's currently at 03.19.00) and change the version number to 03.19.01. Check out / pull dev, then go to your branch and run git merge dev.

Also, a quick reminder on one optional point to consider while you are doing this:

Just one question. WhiteDwarfs.h now includes the following comment next to the declaration of DetermineAccretionRegime():
"Can also change flags related to SN events"
Does this refer only to m_HeShellDetonation, or some other flags as well? Might be useful to spell this out.

@jeffriley , I think you'll need to approve this as well since you requested changes earlier.

[jeffriley]

Thanks @nrsegovia and @SimonStevenson

With the constants added to constants.h, under the comment

// Constants for WD evolution

Can you just ensure that it is clear where the constants come from? I know the constants are named, but is it always clear what paper the constants came from? Maybe, if necessary, just put a comment prior to each group of constants indicating the paper they were taken from?

Modulo that, I approve. I'll mark the PR as approved, and if you can add papers to constants.h (if necessary) before merging, that'd be great.

Also see @ilyamandel's most recent comments. Wrt merging in dev and changing the version number, there is also a conflict with what's new since it was modified in v03.19.00

[SimonStevenson]

Thanks for signing off Jeff.

I have merged in dev and fixed the conflicts in the changelog and What's New page. I also added the references for the WD constants. Does that all look ok @nrsegovia ?

Also can you respond to Ilya's query above:

Just one question. WhiteDwarfs.h now includes the following comment next to the declaration of DetermineAccretionRegime():
"Can also change flags related to SN events"
Does this refer only to m_HeShellDetonation, or some other flags as well? Might be useful to spell this out.

Once this is addressed I think this is ready to go.

[ilyamandel]

Thank you, @SimonStevenson , looks perfect -- and ready for a merge once @nrsegovia gives it a final check.

[nrsegovia]

Apologies for the delay. @SimonStevenson the references in constants.h look good to me. About @ilyamandel's comment: perhaps the statement can be improved? The two flags that can be modified are related to off-center ignition (CO WD evolution into ONeWD) and He-shell detonation, so we could spell it out. But looking at the related code, I wonder if He-shell detonation is properly implemented: it seems that ResolveHeSD() is never called, unlike ResolveSNIa(). Can someone confirm this before I change the code?

[SimonStevenson]

Indeed it does seem that ResolveHeSD is never called.

[ilyamandel]

@nrsegovia :

m_HeShellDetonation is used in COWD::IsSupernova() while m_OffCenterIgnition is used in COWD::EvolveOnPhase() [see COWD.h for both]. But my point here was just that we could explicitly list the flags that are changed. This function doesn't need to describe how they may be used elsewhere.

The fact that we have a function that's not used at all -- ResolveHeSD() -- is a separate issue that I hadn't appreciated. If we don't use it, we should remove it.

[ilyamandel]

I updated the comment. Will merge now, since ResolveHeSD() is a separate issue. @nrsegovia , would appreciate your thoughts on whether that can be safely removed.