Follow-on issues from GW radiation update

[ilyamandel]

There are a few follow-on issues from #1080 (not problems with that PR itself, but things that became apparent while working on it) that I am listing below so we don't forget them:

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. 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. Otherwise, we end up having to use hacks such as setting the semi-major axis to small positive numbers as the binary approaches GW-driven merger.

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.

4. Now that we have a new ChooseTimestep() function in BaseBinaryStar to account for GWs, we should use the same functionality for tides.

[veome22]

@ilyamandel I would be happy to add the timestep code for tides. Quick question-- what would you think about defining global variables such m_DOmega1Dt_tidal, m_DOmega2Dt_tidal, etc. as in the GW code? I've been avoiding it, but it would reduce the number of computations per timestep.

[ilyamandel]

Thank you, @veome22 !
To be honest, I can't recall why global variables were really necessary for either tides or GWs -- why can't the functions that compute them just return them to the caller in a tuple? -- but if you think they are needed, go for it.

[veome22]

Hi @ilyamandel. My current internal implementation of tidal timesteps calls the relevant functions, but the evolution can become painfully slow for the strongest tides. I expect that calculating everything only once and then storing as global variables would increase overhead, but hopefully cut down the run time. I can test both implementations and let you know the time difference, if any.

[ilyamandel]

The fourth item has been addressed, so un-assigning @veome22 since the first three are in my court. :)

[ilyamandel]

@nrsegovia writes the following:

I think there are also some "not-yet-properly-handled" events when both the option to let double WD systems evolve and the relatively new GW radiation prescription are enabled. I have noticed that sometimes a WD would overflow its Roche Lobe (OK), but then lose almost its entire mass. The companion does not experience any mass changes, though the orbit is properly updated. I think this "might" be okay as the companion could be accreting and experiencing flashes ("zero net accretion", no mass changes), but I am not sure about how to deal with the donor... it feels wrong to continue calling it a WD if it loses most of its mass. Probably related to this issue is the fact that such mass transfer takes a while, and results in really long detailed evolutionary histories (~600k steps with the default time step configuration). It also slows down the overall run time of a given population.

First of all, we should check that mass transfer from a WD can actually remain dynamically stable, since donor WDs ought to re-expand on mass loss (maybe we don't have that properly implemented?). Secondly, if the mass transfer is truly stable (the donor can stay within its Roche lobe, the accretor can process the mass), then I would expect the mass losing donor to remain a WD. Do you disagree, @nrsegovia ? Is this something you are willing to look into?

[nrsegovia]

@nrsegovia writes the following:

I think there are also some "not-yet-properly-handled" events when both the option to let double WD systems evolve and the relatively new GW radiation prescription are enabled. I have noticed that sometimes a WD would overflow its Roche Lobe (OK), but then lose almost its entire mass. The companion does not experience any mass changes, though the orbit is properly updated. I think this "might" be okay as the companion could be accreting and experiencing flashes ("zero net accretion", no mass changes), but I am not sure about how to deal with the donor... it feels wrong to continue calling it a WD if it loses most of its mass. Probably related to this issue is the fact that such mass transfer takes a while, and results in really long detailed evolutionary histories (~600k steps with the default time step configuration). It also slows down the overall run time of a given population.

First of all, we should check that mass transfer from a WD can actually remain dynamically stable, since donor WDs ought to re-expand on mass loss (maybe we don't have that properly implemented?). Secondly, if the mass transfer is truly stable (the donor can stay within its Roche lobe, the accretor can process the mass), then I would expect the mass losing donor to remain a WD. Do you disagree, @nrsegovia ? Is this something you are willing to look into?

For reference, I found this behavior when looking into the following system:

COMPAS --allow-touching-at-birth TRUE  --detailed-output TRUE -n 1 -z 0.011 --initial-mass-1 4.139132698621 --initial-mass-2 1.181877457321 --semi-major-axis 0.194801814854 --emit-gravitational-radiation True --evolve-double-white-dwarfs True 

~594,000 steps for the last ~7 Myr of evolution, seems a bit odd. Though I admit that there are other parts of the evolutionary history that are odd, such as the secondary producing a ~0.2 HeWD (I guess it overflowed its Roche Lobe early on the FGB stage?) and then losing most of its mass in a second overflow episode.

Please ignore the (1) on the y-axis label.

Around the time this happens, there is also a "bounce" in the orbital separation, I do not know if this is expected. The separation is really small, too (note the log in the following plot):

What I meant by not being sure about the donor still being a WD is that I'm skeptical about the final ~0.014Msun mass of the donor. I can try looking some more into this issue, but I'm also okay with someone else looking into it full time.

[ilyamandel]

Hi @nrsegovia ,
I am happy to look at the code, but you know more about white dwarfs than I do. The white dwarf radius in COMPAS is given by Eq. (91) of Hurley+ (2000), which I copied below. I have checked that it's correctly implemented in WhiteDwarfs::CalculateRadiusOnPhase_Static(). At the moment, we allow white dwarfs to donate mass on their thermal timescale. We check if there's a solution where the white dwarf would fit inside its Roche lobe after mass loss. When the WD's mass becomes very low, the mass ratio is very extreme, so mass loss widens the binary enough to in principle allow some tiny remnant of the WD to fit into its Roche lobe. Is that physical? If not, should we put in some threshold? E.g., if more than half [or some other fraction] of the WD's mass would need to be lost before it can fit into the Roche lobe, this is a merger?

[nrsegovia]

The white dwarf radius in COMPAS is given by Eq. (91) of Hurley+ (2000), which I copied below.

@ilyamandel I believe that the continuous RLOF might be unphysical, particularly due to the low mass of the remnant. Eq. 91 in Hurley+ (2000) seems to diverge as it approaches to 0, so I would think that this would be understood as a continuous state of RLOF by the code. An alternative treatment is shown in Marsh+ (2004), which can be traced back to Eggleton (1986, though it is referenced as priv. comm.). See the M vs R plot below:

From Marsh+ (2004):

The advantage of this relation is that it provides us with a single relation that matches Nauenberg’s (1972) high-mass relation but also allows for the change to a constant density configuration at low masses (Zapolsky & Salpeter 1969). Eggleton’s relation applies for the whole range 0 < M < M_Ch...

We did not have any problems before because either the evolution of WDs was halted at double WD formation or, when it was not, we did not have GWR to allow further orbital evolution (and mass transfer). This resulted in WDs being assigned constant masses above ~0.3Msun since birth, with radii that would not be extremely large.

E.g., if more than half [or some other fraction] of the WD's mass would need to be lost before it can fit into the Roche lobe, this is a merger?

Marsh+ (2004) also discuss mass transfer stability (though also see Piersanti+2015 and Chen+2022), so I would like to take some more time to look into this.

Anyway, I do not think it is related to issues with GWR anymore, so I will open another issue while I work on these improvements 👍.

[ilyamandel]

Thank you for investigating, @nrsegovia ! I agree, there is no "problem" in the sense that the code is doing what we asked it to do -- it's just that we may need to change what we ask in terms of the WD mass-radius relation and/or the stability of mass transfer from WD donors. I look forward to hearing your thoughts on this when you are ready (no rush!).

[ilyamandel]

@jeffriley -- you may be interested in some of the points raised in this issue in view of your current work on one-step-at-a-time