Improvements to WD evolution

[nrsegovia]

Now that GWR has been implemented, further evolution of double WD systems and the overall details of WDs' response to mass transfer events should be properly assessed and improved wherever possible. I think that this should also be good for the next methods paper. I'm opening this issue to allow discussion related to WDs, though some details have already been mentioned in #1199 and #1343 .

Known problems

• Radial response to mass changes: Currently WDs' radii depend on mass, following Eqn. 91 of Hurley+ 2000:

This evidently diverges when mass gets close to 0. An improvement would be to use Eggleton's equation as quoted in Marsh+ 2004. Below you can find the equation and a comparison between the two approaches:

I would argue that using Eggleton's equation instead of Hurley+ 2000 is the way to go, but I'm open to adding it as an alternative prescription if needed.

• AIC and SNe events: Issue WD SNe are not recorded in the Supernovae output files #1350 was created some time ago, but @SimonStevenson has also commented that it has been quite difficult to find a system that undergoes AIC (meaning that we might have problems beyond recording events). While this could be expected due to the specific requirements for an AIC (ONeWD accreting until it goes over Chandrasekhar mass), these events need further testing.
• Which also brings us to mass transfer stability. The behavior of all prescriptions must be tested, though I have already noticed that the qCrits from the HURLEY_HJELLMING_WEBBINK were mistakenly set to 0 (my bad), and @reinhold-willcox has commented that the Ge et al. prescription could be tweaked a bit for WDs.
• And finally, @SimonStevenson also noted that a system that shows a pre-AIC configuration (or similar) does not show accretion onto the primary, even though the secondary does transfer mass. I wonder if this is related to mass transfer handling in stars with no envelope, or am I looking at the wrong code block? For these systems we are interested in stable mass transfer.

I think I have found a ONe-WD - CO WD binary that starts mass transferring (it seems to be stable too), but the ONeWD doesn't seem to accrete much (any?) mass during the mass transfer (and doesn't undergo AIC).
./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
The secondary WD (the donor) does end up losing most of its mass and ends up with a very low mass (0.0097 Msun), so the mass-radius relation may well be an issue here as well.

Please flag any other issues related to WDs here.

[SimonStevenson]

Thanks for starting this issue Nicolas! If you have a working implementation of the Eggleton WD mass-radius relation, I'm happy to help test/review that. As you say, I think it probably makes sense to add it as an alternate option.

I'll look into what happens in that mass transfer code block and report back.

[SimonStevenson]

The immediate culprit for the ONeWD not accreting seems to be that maximumAccretionRate and betaThermal are both 0, leading to FractionAccreted = 0 and massGainAccretor = 0. It seems that we do not have CalculateMassAcceptanceRate implemented for ONeWDs (unlike for HeWD and COWDs, which each have their own implementations). ONeWD does not inherit from COWD, but from WhiteDwarfs, which inherits from Remnants. Ultimately, I'm not sure which version of CalculateMassAcceptanceRate it ends up using.

I've added a version of ONe::CalculateMassAcceptanceRate to ONeWD.cpp, but am not completely sure what assumptions we should be making there. @nrsegovia Does it make sense to use the same as for COWD (following Claeys et al), or something different?

[nrsegovia]

Thanks for starting this issue Nicolas! If you have a working implementation of the Eggleton WD mass-radius relation, I'm happy to help test/review that. As you say, I think it probably makes sense to add it as an alternate option.

I came up with the following:

double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) {
    //From Eggleton 1986, as quoted in Verbunt & Rappaport 1988 and Marsh et al. 2004
    double MP = 5.7e-4; // Constant 
    // sanity check for mass - just return 0.0 if mass <= 0
    if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0;
    
    double MCH_Mass_one_third  = std::cbrt(MCH / p_Mass); 
    double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third;
    double MP_Mass = MP / p_Mass;
    double MP_Mass_two_thirds = MP_Mass / std::cbrt(MP / p_Mass); 
    double First_Factor = std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds));
    double Pre_Second_Factor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass;
    double Second_Factor = std::cbrt(Pre_Second_Factor) / Pre_Second_Factor;

    return 0.0114 * First_Factor * Second_Factor;
}

It uses the same CalculateRadiusOnPhase_Static as my intention was to supersede the previous implementation. Also, note that the plot I had in my initial comment has been edited. I had incorrectly assigned an exponent.

The immediate culprit for the ONeWD not accreting seems to be that maximumAccretionRate and betaThermal are both 0, leading to FractionAccreted = 0 and massGainAccretor = 0. It seems that we do not have CalculateMassAcceptanceRate implemented for ONeWDs (unlike for HeWD and COWDs, which each have their own implementations). ONeWD does not inherit from COWD, but from WhiteDwarfs, which inherits from Remnants. Ultimately, I'm not sure which version of CalculateMassAcceptanceRate it ends up using.

I've added a version of ONe::CalculateMassAcceptanceRate to ONeWD.cpp, but am not completely sure what assumptions we should be making there. @nrsegovia Does it make sense to use the same as for COWD (following Claeys et al), or something different?

I think that here we have two options, either look for studies specifically centered around ONe WDs' response to mass accumulation, or default to the same CO WD treatment. The former includes Guo+ 2021, but I did not find any other recent+relevant studies (and this one only deals with He-rich material). Defaulting to CO WDs' behavior seems reasonable, and it is what was done in StarTrack (section 5.7.3 of Belczynski 2008).

As a side note, and with #209 in mind, it might also be a good idea to eventually consider merger products as in StarTrack (section 5.7):

If dynamical instability is encountered for a binary with two WDs, we assume that a merger occurs. Based on the results of Saio & Nomoto (1998) mergers of a ONe WD with any type of WD companion and two CO WDs lead to either AIC and NS formation (if the total merger mass Mmerger is above MECS = 1.38 Msun ) or the formation of the single ONe WD (with new mass equal to Mmerger). For mergers of CO WD and He WD, we assume a Type Ia SN explosion: either sub-Chandrasekhar (Mmerger < 1.44 Msun ) (see Woosley et al. 1986; Woosley & Weaver 1994) or Chandrasekhar mass SN Ia (Mmerger > 1.44 Msun).

[SimonStevenson]

I've tried using the COWD version of CalculateMassAcceptanceRate for ONeWD and FractionAccreted becomes non-zero, but very tiny (order 1e-31). So the ONeWD still accretes virtually no mass. I'm not really sure why it is so small.

There is also Brooks et al. 2017 https://iopscience.iop.org/article/10.3847/1538-4357/aa79a6 which is also for He accretion on to ONeWDs

[nrsegovia]

I think I might have spotted a bug. Can you check what happens if you change double logMassTransferRate = log10(p_MassTransferRate); to double logMassTransferRate = log10(p_MassTransferRate / MYR_TO_YEAR); in both WhiteDwarfs::CalculateEtaHe and WhiteDwarfs::CalculateEtaH? I think that, at the very least, the EtaHe regime comparisons should be based on Msun/yr rates, while the mass transfer rates are handled as Msun/Myr in COMPAS (please correct me if I'm wrong). This was correctly implemented in COWD::DetermineAccretionRegime but not in the calculation of the eta values. Whoops.

Another option would be to artificially set WhiteDwarfs::CalculateEtaHe to always return 1.0, because in the end this sets the accretion rate. If the donor is not helium rich, then you would also have to set WhiteDwarfs::CalculateEtaH to test whether the problem stems from here.

[SimonStevenson]

Well spotted @nrsegovia ! Making that change caused the ONeWD to collapse to a NS.

The SN is correctly recorded in BSE_Supernovae as SN_Type(SN) 32, which corresponds to AIC. The mass of the NS ends up being 1.26 (the M_ECSN), which is less than the progenitor ONeWD mass.

The system still has this 'separation bounce' later on and experiences another phase of MT as a NS-COWD binary. I'm not sure if that is physical, or should just be a merger.

So changes for the PR:

• Add ONe::CalculateMassAcceptanceRate . For now, just using same as for COWD.
• Fix units of logMassTransferRate in WhiteDwarfs::CalculateEtaHe and WhiteDwarfs::CalculateEtaH
• Update white dwarf mass-radius relation (WhiteDwarfs::CalculateRadiusOnPhase_Static)

[SimonStevenson]

I can put together the PR during the code call later today

[nrsegovia]

That sounds good to me, it should take care of most of the problems stated on the first post. I would like to make some (minimal) additional changes to the documentation and the qCrits.

About the separation bounce, I was skeptical at first, too. However, it seems to be something expected in mass transferring compact binaries such as cataclysmic variables. I'm slightly concerned about the mass ratio being a bit "off":

Mass transfer between two white dwarfs has generally been taken to be stable if the mass ratio q = M_2/M_1 , where M_1 is the mass of the accretor, is smaller than about 2/3.

This is from Marsh+ 2004, while Hurley+ 2002 quote a slightly lower number, 0.628. So perhaps it would be a merger? I'm not sure about how to evaluate this for the default zeta prescription.

[ilyamandel]

I chatted a bit with Alexey Bobrick today about WDs as donors.

Alexey reminded me that WDs do have masses of 0.0few Msun in ultra-compact X-ray binaries. So, indeed, our general model of WDs getting rapidly stripped until the mass ratio become very extreme, at which point the system can become a relatively long-lived system (a bit like an LMXB, but with driving from gravitational waves rather than nuclear evolution), is sensible. Now, we currently assume that this happens on the thermal timescale of the donor -- which for remnants is set to the dynamical timescale, which would only be tens of seconds for a WD. That may be the response timescale of the WD to stripping, but detailed calculations show that the actual timescale of the rapid phase of mass transfer is rather longer, hours-days-weeks. However, since that's still much faster than any reasonable timescales of the accretor, this shouldn't make any practical difference.