Convergence issues starting from the HeMS

[reinhold-willcox]

Describe the bug

The outcomes of stars that spend some time as HeMS stars is showing some convergence issues. I find that the He and CO core masses vary by a lot (in some cases by a factor of ~2), depending on how many timesteps are resolved during the HeMS phase. This may be occuring if the wind mass loss rate grows substantially along the HeMS, such that the integrated quantity of mass loss is much larger if the Delta T is reduced, compared to the case where there is only one timestep covering the entire HeMS.

In the plots below, I show a small sample of very wide (effectively single) binaries, with a companion chosen for discrete mass ratio values between 0.1 and 1. I find that when the mass ratios are more equal (q > 0.5), the final core masses diverge. For lower mass companions, their nuclear lifetimes are longer, so they are not the limiting factor in the calculation of timestep duration. For higher mass companions, the opposite is true, and they reduce the timestep duration more than that required for the primary. As a result, the primary spends more time on the HeMS, resulting in artificially greater mass loss from winds.

When I reduce the timestep by a factor of 10 across the board, the outcomes converge much better for different values of the mass ratio, but they converge to final core masses that are even significantly lower than in any of the cases with normal timestepping. This is an unphysical artefact of undersampling in this phase. I find that this starts to occur at M_ZAMS ~ 40 Msun, meaning most of our conclusions related to core masses and black holes from progenitors above this threshold are potentially affected.

Initial M_ZAMS to final M_CO relation, for a population of effectively single stars in wide binaries at Zsol, colored by mass ratio.

Late stage evolution of identical 140 Msun stars at Zsol. Colors are same as above, you can see that the M_He diverges before CC for the standard timesteps. There are also no timesteps where M_CO is non-zero before collapsing into a BH, but the BH masses are set here by the He core mass so they too are divergent.

Final He and CO core masses (for the same 140 Msun stars as above) as a function of q, for standard and reduced timesteps.

Label the issue
urgency_moderate - This is a moderately urgent issue
severity_major - This is a severe bug

Expected behavior
Slight differences in the final masses can be expected any time you discretize a continuous process, but differences by a factor ~2 suggest we generally need much smaller timesteps along the HeMS.

Versioning:

• COMPAS v03.07.00

[jeffriley]

Also see #24

[SimonStevenson]

This is rather concerning! Your set up is a little artificial though (not our typical use case). Is this using the default values for --mass-change-fraction and --radial-change-fraction, or --mass-change-fraction 0.005 --radial-change-fraction 0.005 or some other values (it would be good to include an example command to reproduce this)? Or did you use --timestep-multiplier? We know that the defaults do result in some convergence issues (though I've never seen this one). Do you know if this issue similarly manifests in a typical close interacting binary (e.g., GW progenitor)?

[reinhold-willcox]

@SimonStevenson Sorry for the delay, I missed your message.

These variations used only --timestep-multiplier, not the --X-change-fraction's. I can try varying them as well though it may take a while before I have some bandwidth for it.

These changes should affect close binaries as well, since it's really just a timestep integration issue. The problem is mediated, to some extent, for closer binaries, since any RLOF will in a sense add "extra" timesteps, which will help the integration to be more accurate. I think the solution is just to require more timesteps along the HeMS (even though this will increase run time per binary)