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.

[ilyamandel]

I think it would be easier to understand what goes wrong by looking at the evolution of a single binary rather than the end-point of a whole population.

In principle, even when we are transferring a lot of mass in one step, we do split up the transferred mass into lots of small pieces in BaseBinaryStar::CalculateMassTransferOrbit(), and the code there looks sensible -- but it would be really good to understand if we are doing something wrong!

[ilyamandel]

@SimonStevenson, @yuzhesong : OK, I checked the code for mass transfer, and I don't see any evidence of bugs -- it seems to be doing what was intended. So, the question is whether the intent was correct. :)

There are two places where one might quibble with the code assumptions.

1. We assume a fixed mass transfer rate for the entire episode as computed in a given time step, which is pre-computed at the start of the timestep. That means both donor and accretor properties are considered only at the start of the timestep, even if we donate, say, the entire envelope of an evolved star on a single timestep. Note that we correctly integrate the orbital evolution by taking small substeps, but the accretion efficiency beta is determined once and for all. Given that we are using the thermal timescale for evolved stars, I expect this might change the donor's nominal thermal timescale (and hence, accretion efficiency) by a factor of ~2 [R is roughly constant at the RL size, L is constant, but M can ~halve]. However, given that the accretor is changing even more drastically and our treatment is an over-simplification regardless (see, e.g., Lau+ 2024 on Hamstars, which we should eventually move to using), I don't think this is the key issue.

2. A larger issue for mass transfer from MS stars onto compact objects is that we are using the donor's thermal timescale to estimate mass transfer efficiency, comparing it with the accretor's (Eddington-limited, in this case) accretion timescale. That generally leads to extremely non-conservative mass transfer. However, the nuclear can be much longer than the thermal timescale for MS stars, particularly low-mass MS stars, i.e., closer to the accretion timescale, and hence may allow for much slower mass transfer. I suspect this is the also responsible for #889 .

I propose that we leave choice 1 as is, but work on choice 2.

[ilyamandel]

OK, I have a plane for how to allow for nuclear timescale mass transfer without taking too much of a hit on computational cost, will submit a PR either today or Monday (offline for the weekend) and ask you and Jeff to review.

[ilyamandel]

Should be at least partially addressed in #1123 , which allows MS stars to donate mass on their nuclear timescale (see above). @SimonStevenson , @yuzhesong , let me know how things go once that's pulled in (not closing for now since I haven't tested consequences specifically for NS accretion).

[ilyamandel]

Sorry, I spoke prematurely, and there were still significant problems with case A mass transfer following the previous PR. Please try #1132 instead (once it's reviewed and pulled in).

[ilyamandel]

Somehow, the relevant change didn't seem to make it into #1132 (mysteries of git?), but should be there in #1133, so please use that instead, once approved.

[ilyamandel]

@yuzhesong , @SimonStevenson -- have you been able to check if this issue is now resolved? So far, I seem to find that the MT is behaving much better than previously in all of the tests I've considered, including accretion onto compact objects.

[yuzhesong]

@ilyamandel your fix in #1133 does indeed let us have long case A mass transfer, but it seems like the mass transfer rate is too low. I made a comparison between your fix and including your fix and Simon's treatment in the first figure below.

This seems to indicate that NS won't accrete enough mass to let it reach to the typical MSP magnetic field strength at around 10e8 Gauss, shown below.

And Figs. 2a, 2b, 2c & 2d in this Tauris et al 1999 paper (https://ui.adsabs.harvard.edu/abs/1999A%26A...350..928T/abstract) also seems to indicate a higher mass transfer rate.

We will run some MESA models for some more insights.

[ilyamandel]

@yuzhesong : The nuclear timescale of a 0.9 Msol star is more than 10^10 years. That means I don't expect it to transfer mass faster than 0.9 Msol / 10^10 years, or 10^{-10} Msol per year = 10^{-4} Msol / Myr.

The Mdot plot you made shows that we currently estimate it as being a bit over 10^{-5} Msol / Myr, which I consider sensible (we are using the radial expansion timescale, which is lower than 10^10 years except late on the MS).

I am not sure how @SimonStevenson arrived at 5e-3 Msol/Myr, but that seems to me to be rather high [if this is taken from BSE, I think BSE rates are artificial because they depend on the amount by which the star exceeds the Roche lobe, which just depends on the size of the timestep].

Do you have a physically motivated alternative prescription for nuclear timescale mass transfer? We can either replace the current one, or create an alternative option, which would certainly be in the spirit of rapid pop synth. Alternatively, we can also parametrise our ignorance by including a new option as a multiplier to the existing nuclear timescale mass transfer rate.

[ilyamandel]

Aside: correction to MT in #1147, please use that for further testing.

[yuzhesong]

@ilyamandel noted, will test and if everything works OK, then we'll merge the MSP branch into dev.

[SimonStevenson]

Just wanted to follow up on this issue. I have tried redoing the comparison to Tauris & Savonije. things look a bit better, but still do not match what they get.

It seems that the mass transfer now happens all in one timestep, rather than over a a prolonged period of time as above

(I'm also not sure why there are two clusters of lines)

As another data point, we got a plot from Jan showing what BPASS gets, and they also reproduce the Tauris & Savonije relation.

As a reminder, I am generating the range of post-CE binaries by adjusting the common envelope efficiency with something like:

./COMPAS --common-envelope-alpha s[1.0,1.1,1.2,1.3,1.5,1.7,2.0,2.3,2.5,2.7,3.0,4.0,5.0,7.0,10.0] -n 1 --evolve-pulsars TRUE --detailed-output TRUE --initial-mass-1 10.0 --initial-mass-2 1.0 --orbital-period 1000.0 --chemically-homogeneous-evolution-mode NONE --maximum-evolution-time 13700 --radial-change-fraction 0.01 --mass-change-fraction 0.01 --kick-magnitude-distribution ZERO --remnant-mass-prescription Fryer2012

[ilyamandel]

@SimonStevenson -- I am not sure from whether it's our MT that's somehow different from Tauris & Savonije, or our method for calculating the final core/WD mass that's somehow different. I've tested the former quite a bit, so I am, for the moment, more suspicious of the latter.

Figure 5 of Tauris & Savonije that you posted some messages back seems to suggest that the core mass of the donor increases even as it loses mass -- it moves up and to the right, toward larger orbital periods and higher core masses. Is that what happens in COMPAS? And, if so, do the core masses match, or are the COMPAS purple dots on your last population-level plot just shifted leftwards from where they should be?

Can you identify a couple of COMPAS tracks to directly compare to figure 5 of Tauris & Savonije?

[SimonStevenson]

I tried making a plot showing the evolution of an individual system, but it looks kind of hilariously bad:

(M = helium core mass here)

The post-CE/SN orbital period is around 2d, and mass transfer starts when the 1 Msun companion has a core mass of around 0.16 Msun. This isn't too far from the closest model in Tauris & Savonije. In their model, the core mass continues to grow to 0.24 Msun, which doesn't happen in the COMPAS model. So that would move the COMPAS point to the right a bit.

The other issue seems to be that in the Tauris & Savonije model the mass transfer is almost completely conservative (see the bottom middle panel of Fig 2a), so the orbit only grows to ~10 d, rather than 100d in the non-conservative COMPAS model (which only accretes 0.025 Msun out of ~0.8 Msun transferred).

So it might be that the COMPAS model needs to move both right and down.

I think part of the problem may be that the helium core grows over a significant time (~1 Gyr for the 1 Msun model), and our mass transfer is treated as an instantaneous process using whatever the core mass was when the binary fills its Roche lobe. For example, I made this plot for a grid of single stars (0.5--1.5 Msun) showing the growth of the helium core mass with time. Only the stars with initial masses > 1 Msun develop cores within the maximum evolution time.

[I did some more reading on the Porb-MWD relationship. It seems that this Rappaport et al. 1995 paper is an earlier reference for the same relation: https://adsabs.harvard.edu/full/1995MNRAS.273..731R -- see also these Nature papers on the formation PSR B: https://ui.adsabs.harvard.edu/abs/1983Natur.304..419J/abstract Paczynski 1983 https://ui.adsabs.harvard.edu/abs/1983Natur.304..421P/abstract Savonije 1983 https://ui.adsabs.harvard.edu/abs/1983Natur.304..422S/abstract]

[ilyamandel]

Thank you, @SimonStevenson . It sounds like Tauris & Savonije think this should be nuclear timescale MT while COMPAS thinks it should be thermal timescale. Can you please share the parameters of this binary (i.e., the command line)? I'd like to explore it in a bit more detail.

[SimonStevenson]

That sounds like it could be the difference, and I imagine that can be fairly sensitive to exactly when the companion fills its Roche lobe, so a small difference in orbital period could lead to a large difference in evolution. It's hard to get an exact match, which is why I'm using the common envelope alpha parameter to adjust the post-CE separation.

You can reproduce this system with

./COMPAS --common-envelope-alpha 1.2 -n 1 --evolve-pulsars TRUE --detailed-output TRUE --initial-mass-1 10.0 --initial-mass-2 1.0 --orbital-period 1000.0 --chemically-homogeneous-evolution-mode NONE --maximum-evolution-time 13700 --radial-change-fraction 0.01 --mass-change-fraction 0.01 --kick-magnitude-distribution ZERO --remnant-mass-prescription Fryer2012 

[ilyamandel]

Hi @SimonStevenson , Sorry it took me so long to return to this. I find that in your example (running with COMPAS v03.08.04), the first and only mass transfer from the secondary to the primary happens when the secondary is an FGB star. Now, at the moment, we only have nuclear timescale MT enabled for MS donors -- we assume that evolved donors always lose their entire envelope on MT.
Is that a reasonable assumption? Do Tauris and Savonije or other authors allow for nuclear timescale MT from MS donors? If so, how do they evaluate this nuclear timescale?

[ilyamandel]

I recommend testing with #1287 , which makes a number of further improvements to MT treatment.