Issue with mass transfer onto neutron star

[yuzhesong]

Describe the bug
When attempting to fix the issues with neutron star recycling in the current version, @SimonStevenson and I have discovered that in some situations, the NS would accrete 0.25 Msol in one timestep of 1000 years, which is unphysical, since the accretion rate onto a neutron star should be Eddington limited.

@SimonStevenson and I have narrowed down the possible issue in the function
void BaseBinaryStar::CalculateMassTransfer(),
where m_FractionAccreted is calculated using CalculateMassAcceptanceRate().

It seems like that the code is calculating the mass transfer rate in this timestep described above as thermal mass loss rate, while it is much shorter than thermal time scale.

Label the issue

urgency_moderate - This is a moderately urgent issue. While this creates some issues, these cases are relatively rare and don't seem to pose problems to the majority of simulations.

To Reproduce
I have attached a yaml file (in the format of txt as GitHub doesn't support uploading yaml file) that should be able to reproduce this case
config_accretion_bug.txt

Versioning (please complete the following information):

• COMPAS v02.41.05

[SimonStevenson]

I think I have a solution for this. I think the issue is that for these low-mass companion stars, the duration of the mass transfer should be long as the stars have long thermal timescales tau_th (tens of Myrs), which are much longer than the timestep dt we were taking during the mass transfer (in this case, the initial tilmestep is as small as dt = 0.001 Myr, since we are trying to resolve the early pulsar evolution).

To attempt to solve this, I have added a clause that if the star's thermal timescale is longer than the current timestep, we only accrete a fraction (dt / tau_th) of the mass that the donor star wants to donate.

This seems to have the desired effect of avoiding a large amount of mass being accreted in one timestep.

When I first did this, the system ended up merging after a few time steps. I managed to track this down to the tolerance ROOT_ABS_TOLERANCE used in MassLossToFitInsideRocheLobe, which by default is 1E-10. I noticed that the solution that was found was good to a few 1E-10. Relaxing the tolerance to 1E-9 solved this issue.

[ilyamandel]

Hi @SimonStevenson ,
I think we speculated that this might be the issue during one of the code calls. However, it seems that the NS is actually doing the right thing: it is using the correct mass transfer rate (slow thermal timescale, not delta M / delta t (where the latter is the time step size) to decide that it can, in fact, accrete a sufficient fraction of mass while remaining sub-Eddington. Is that correct?

[SimonStevenson]

@ilyamandel that's right, with this fix, the system now evolves correctly. This will be part of @yuzhesong 's PR to fix MSP formation. Unfortunately, it seems that this issue (or one similar to it) persists even with this fix.

Here the different coloured lines refer to different values of alpha common envelope, which I am using as a proxy to manually adjust the post-common envelope orbital separation/orbital period. I looked at the relation between the final He WD mass and the final orbital period (at 14 Gyr) and compared this to the relation from Tauris & Savonije. It seems the relation is very different in COMPAS.

Compare this to this figure from Tauris & Savonije:

[ilyamandel]

OK, I think there are a couple of things here, @SimonStevenson :

• I didn't mean that we need a "fix" for the MT, I actually meant that we seem to be doing the right thing (in terms of how conservative MT should be overall), but maybe you don't want too much to be transferred in one time step because you are concerned that there are other driving forces at play? The reason we tend to transfer a lot of mass at once is that we don't really know how to describe stars that are actively mass-transferring in COMPAS..
• Are there other drivers besides mass transfer in Tauris & Savonije? E.g., magnetic braking or gravitational waves? Both are on the to-do list (Adam has a GW implementation, just waiting for some fixes), but not yet in. I ask because their binary continues hardening as the low-mass donor is mass -- you need something else to bring the binary in, mass transfer should widen it as it does in COMPAS models.

[SimonStevenson]

1. Well there's a big difference between a system that transfers 0.25 Msol in 1000 years and one that transfers a similar amount of mass in ~600 Myr. Not least of which is that it was causing problems for our pulsar recycling code, but any attempt to post process these systems for LMXB/MSP studies will result in the wrong answers for lifetimes/birthtimes etc. Switching between main sequence models here should be ~fine, I think.

2. Yes they do include those, and for close binaries (with Porb < 1 day) those effects can certainly become important, but for the wide binaries with Porb >> 1 day studied in Tauris & Savonije, they should be mostly negligible (and they say that in their paper). So all that should be important here is the nuclear evolution of the low-mass companion star that drives the mass transfer (which should be ok in COMPAS) and the details of mass transfer.

[ilyamandel]

OK, fair point on 1. Then we'll need to think about how to deal with MT more broadly, if we want to split up things like, say, the fast phase of case A MT over multiple time steps...

On 2, I take your point, but then don't understand how their binary shrinks on mass transfer from a low mass WD to a high mass NS. One can analytically confirm that both fully conservative and fully non-conservative (isotropic re-emission from accretor) MT should widen such binaries. Do they assume very large angular momentum carried away by non-accreted mass?

[SimonStevenson]

They make the standard assumption of isotropic re-emission from the accreting neutron star. The binaries do widen. They start with orbits of a few days and end with orbits of 2--1000 days. Sorry if the plot above was confusing, that's not an evolutionary track, but a locus connecting the final evolutionary states for a series of different initial configurations. That can be seen more clearly in this figure:

[SimonStevenson]

I ran a similar experiment with BSE and it seems to get closer to the Tauris & Savonije result, though still not quite the same relation.