Update pulsar accretion treatment in NS.cpp

[yuzhesong]

This PR is created to update to neutron star accretion treatments, to fix MSP formation/NS in mass transfer treatments:
1). Created a new function NS::PulsarAccretion() to calculate the pulsar evolution in stable mass transfer.
2). In UpdateMagneticFieldAndSpin(), splitting stable mass transfer into smaller steps so that no negative spin period is present.
3). Adding a new programing option "NS_ACCRETION_IN_CE" for different treatment of how neutron star would behave when in CE.
4). This is also solution to issue #1002

[jeffriley]

Thanks Robert!

Just quickly for now - I'll review in more detail later:

1. The version number reflects the changes to the code: we use MM.mm.dd, where MM is a major update, mm is a minor update, and dd is a defect repair. You've added new functionality, so this PR is a minor change to functionality, so you should change the version number to 02.43.00
2. You've added a new option, so you need to:
   (a) update the documentation accordingly.
   (b) create a new default yaml file (you can use the "create-yaml-file" option to do that.
3. Your changelog entry notes new option "NS_ACCRETION_IN_CE" - those underscores should be hyphens.

[jeffriley]

@yuzhesong A few more comments:

NS::PulsarAccretion() needs @return in description
NS::PulsarAccretion() parameter kappa should be p_Kappa
NS::PulsarAccretion() has a couple of stray std::cout ... statements - should be removed

NS::UpdateMagneticFieldAndSpin() has a bunch of stray std::cout ... statements - should be removed
NS::UpdateMagneticFieldAndSpin() indentation is out somewhere - lines 635 and 636 have braces in the same column
NS::UpdateMagneticFieldAndSpin() what is the comment at the end of line 599 for? There is no P_f in NS::UpdateMagneticFieldAndSpin()...
NS::UpdateMagneticFieldAndSpin() ditto line 625
NS::UpdateMagneticFieldAndSpin() the comment at the end of line 619 is not useful - remove it

We use camelCase in COMPAS for variable and function names, not snake_case, and we (generally) start local variable names with lowercase, and function names with uppercase. In NS::UpdateMagneticFieldAndSpin(), change:
Accretion_Results to accretionResults
Next_Accretion_Results to nextAccretionResults
divide_timestep_by to divideTimestepBy
do_not_stop_looping to doNotStopLooping
new_timestep_size to newTimestepSize
new_massGain to newMassGain

In NS::UpdateMagneticFieldAndSpin(), this block of code doesn't need to be inside the for loop:

if (n==1) {
    last_B = std::get<0>(Accretion_Results);
    last_f = std::get<1>(Accretion_Results);
    last_fdot = std::get<2>(Accretion_Results); 
    last_am = std::get<3>(Accretion_Results);
}

What is the motivation for setting Next_Accretion_Results = Accretion_Results ; at line 572, when you break out of the for loop at line 573, and Next_Accretion_Results is not used outside that for loop?

Can I suggest you rewrite the existing while loop something like this:

int intraSteps = 100;
bool done      = false;
while (!done) {
    
    intraSteps *= 2;
    
    double thisTimestepSize = p_Stepsize / intraSteps;
    double thisMassGain     = p_MassGainPerTimeStep * 1000.0 / intraSteps; 
    
    accretionResults = PulsarAccretion(m_PulsarDetails.magneticField, m_PulsarDetails.spinFrequency, 
    m_AngularMomentum_CGS, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);

    B    = std::get<0>(accretionResults);
    f    = std::get<1>(accretionResults);
    fdot = std::get<2>(accretionResults); 
    am   = std::get<3>(accretionResults);

    int count = 0;
    while (!done) {
        thisAccretionResults = PulsarAccretion(B, f, am, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);
        f  = std::get<1>(thisAccretionResults);
        
        if (utils::Compare(f, 0.0) < 0) break;

        B    = std::get<0>(thisAccretionResults); 
        fdot = std::get<2>(thisAccretionResults); 
        am   = std::get<3>(thisAccretionResults); 
    
        if (++count >= intraSteps) done = true;        
    }
}

I'll stop here for now and come back to it later.

EDIT:

I forgot to ask - is fdot used anywhere in the code (or going to be)? If not, you don't need to extract it from accretionResults.

And can we put in a guard for the while loop never terminating? Maybe we do the inner loop (or the outer loop - doesn't really matter I guess) a maximum of X times (you pick X)?

[jeffriley]

You can also use std::tie() to receive the results from a function that returns a tuple:

std::tie(B, f, fdot, am) = PulsarAccretion(B, f, am, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);

[yuzhesong]

@jeffriley the last two commits should've addressed all of your comments.

[ilyamandel]

NS.cpp line 388 says " // calculate the spin down rate for isolated neutron stars"
however, the calculation that immediately follows is computing the initial spin, not the spin-down rate -- that follows after line 410, where the same comment is repeated (correctly on that occasion), right?

[ilyamandel]

NS.cpp:

• @return Tuple containing the updated magnetic field strength, spin frequency, angular momentum of neutron star

Actually, you have four returns.

• @param [IN] p_MassGainPerTimeStep Mass loss from the secondary for each iteration (in g)
• @param [IN] p_Kappa Mass loss from the secondary for each iteration (in g)

Repeated comments for different quantities; and I don't think either of the comments is correct or matches the variable name. :)

[ilyamandel]

NS.cpp, line 455: (G * 1000.0)

use G_CGS instead

[ilyamandel]

There seem to be some funny factors of 2 in the Alfven radius. Yours has an extra factor of 2 relative to Eq. 8 of Oslowski, whom you cite; that equation does not match Eq. 1 of Miller, https://www.astro.umd.edu/~miller/teaching/astr498/lecture19.pdf; and I actually don't find a factor of 2 at all when I do it [though I was not very careful]...

[ilyamandel]

double keplerianAngularVelocityAtAlfvenRadius = 4.0 * M_PI * keplerianVelocityAtAlfvenRadius / alfvenRadius;
double velocityDifference                     = keplerianAngularVelocityAtAlfvenRadius - p_SpinFrequency;

// calculate the change in angular momentum due to accretion
// see Equation 12 in arXiv:0805.0059/ Equation 8 in arxiv:1912.02415 
double Jdot        =  p_Epsilon * velocityDifference * alfvenRadius * alfvenRadius * mDot;

That doesn't seem right. First of all, the Keplerian angular velocity is just omega=v/R, no factors of 4 pi. Secondly, if p_Frequency is really the frequency, it should be multiplied by 2 pi to convert it into velocity. Thirdly, the comment for p_Epsilon is misleading. Fourthly, the amount of accreted angular momentum is just p_MassGainPerTimeStep * keplerianVelocityAtAlfvenRadius * alfvenRadius, no need for differences.

[Aside: I now see that you are following Kiel+ in using half the Alfven radius as the magnetosphere coupling radius; fine, but then be consistent with variable names.]

I am stopping here to switch to another project, but suggest going through this and remaining code carefully, making sure equations are correctly transcribed (and correct in the first place), comments and function names are clear, etc.

[ilyamandel]

    // If the current timestep is smaller than the donor's thermal timescale, only a fraction of the mass
    // that needs to be donated has time to be donated
    double donorThermalTimescale = m_Donor->CalculateThermalTimescale();
    if (p_Dt < donorThermalTimescale) {
        massDiffDonor = (p_Dt / donorThermalTimescale) * massDiffDonor;
    }

I worry that this is a very major change with a number of potential unintended consequences. For example, we generally remove the entire envelope mass of an evolved donor in one time step during mass transfer. I am not sure we would correctly treat a donor that effectively pauses its evolution part way through donating its envelope. We would need to do careful regression studies showing the impact of this change on mass transfer in all types of binaries (say, HG onto MS companion, CHeB onto HeMS, etc.), not just onto neutron stars. This should include looking at how we are storing and displaying stellar properties of the donor that is paused part way through thermal-timescale mass transfer in the detailed output files. It's quite a bit of work, but this is a major change that we should not make lightly.

[yuzhesong]

@SimonStevenson and I have been doing some testing on the new MT descriptions (#1125 ) by @ilyamandel and check its impact on neutron star in mass transfer. The figure attached here plots, from top to bottom, using the configuration file attached:

• P-Pdot diagram of a NS in a NS-MS binary system
• Period of NS vs. time
• Pdot of NS vs. time
• Mass of NS vs. time
• Mdot of NS vs. time,
• magnetic field of NS vs. time.

The blue curves (labeled I) in all figure are produced with MT treatments in #1125 , and orange curves (labeled I+S) are produced with implementation in #1125 AND Simon’s treatment on breaking down the nuclear mass transfer of a MS onto NS into smaller timesteps, which is currently implemented in this PR.

Treatments in #1125 still has a one short step of large mass transfer, while the other treatments provide gradual mass transfer. In the end, both treatments transfer the same amount of mass onto the NS, but within different timescale (ass seen in the plot). This means that magnetic field of the NS should decay to the same value (as seen in the plot). But due to the different Mdot (as shown in the second last panel of the plot), the Alfven radius is then different, which then cause the difference in the calculation of P and Pdot of the pulsar at each time step, giving different final results. Overall, it seems like having the mass transfer broken down over time seems a better treatment for pulsar evolution.

Simon and I think that maybe we can use Ilya’s mass transfer treatments in all cases, unless we have evolve-pulsars=True, then we can implement Simon's treatment as a special case for NS in MT?

compasConfig_v2.txt

[ilyamandel]

Hi @yuzhesong ,
Unfortunately, #1125 was not the complete fix I hoped, due to some underlying inconsistencies elsewhere in the code. I hope that #1132 will address the remaining issue.

Just in case, to help me with testing your particular issue, could you give me a couple of COMPAS command lines (i.e., initial conditions for the binaries at ZAMS) that should lead to slow MT onto NSs?

[yuzhesong]

@ilyamandel will wait for #1132 to be merged to check again.

If you use the compacConfig_v2.txt file listed in my reply above, and change it to a yaml file, it should produce a NS-MS binary system in slow case A mass transfer.

[ilyamandel]

The mass of the NS grows from 1.1773 Msol at formation to 1.4322 Msol at the end of the evolution in your example with the latest PR... but then the binary merges, which should not happen. I'll move the PR to draft and investigate ASAP, but probably not before Monday or even Tuesday.

[yuzhesong]

@ilyamandel that is weird and shouldn't have happened. Using the corrections in #1125 and Simon's treatments, this binary system evolves to "allowed time exceeded".

[ilyamandel]

Actually, managed to fix that faster than I thought. You'll be happy to know that we now do, in fact, limit nuclear timescale MT to no more than the MT rate times the timestep duration. However, we only do this for nuclear timescale mass transfer, and we do it for all accretors, so I think this is a more robust solution. I have tested that your NS accretes around ~0.2 Msol (slowly, over many timesteps). I have not tested other features such as the B field strength -- leave that for you, after PR #1132 is reviewed.

[yuzhesong]

@jeffriley @ilyamandel @SimonStevenson

Ilya's new mass transfer fixes work.
Updated the codes for this PR for review to merge to dev.

[jeffriley]

@yuzhesong

"Updated the codes for this PR for review to merge to dev." I don't see any updates - am I missing them? Or are they older changes?

Also, there are conflicts to resolve. And finally - the PR is marked as a draft - do you want to mark it ready for review?

[yuzhesong]

@jeffriley Sorry I didn't manage to push the changes made to PR properly. It should all be resolved now hopefully.

[jeffriley]

All good - thanks @yuzhesong

[ilyamandel]

Commented out code on lines 2145--2151 of BaseBinaryStar.cpp should be removed.

[ilyamandel]

BaseBinaryStar.cpp, line 2177:
// Check for recycled pulsars. Not considering CEE as a way of recycling NSs.
The second part of this comment is no longer correct. The line that immediately follows is a code line that's been commented out -- should just be removed.

@yuzhesong -- just FYI, after creating a PR, I usually go through the "review changes" page myself to see all the places where the code changed, to make sure I haven't left in commented out code, etc.

[ilyamandel]

Options.h:

Line 383 -- what is "neutron-star-in-accretion-in-ce" -- I don't see this mentioned elsewhere?

Line 515: "neutron-star-in-accretion-in-CE" -- this seems to have a different capitalisation than ""neutron-star-in-accretion-in-ce" which appears in Options.cpp

[ilyamandel]

NS.cpp:

line 420: "* The calculations used in this function follows closely Sec. 2.2.1 in arxiv:1912.02415" -- follows->follow

line 465-466:
// see Equation 2 in 1994MNRAS.269..455J / Equation 8 in arxiv:1912.02415
double keplerianVelocityAtMagneticRadius = std::sqrt(4.0 * (G_CGS) * mass_g / alfvenRadius);
I don't know why Eq. 8 is relevant -- maybe you mean Eq. 9?
I don't see why there is a factor of 4 here. If the magnetic radius is 1/2 of the Alfven radius, there should be a factor of 2.
To avoid arithmetical difficulties, just define
double magneticRadius = alfvenRadius / 2.0;
and then use magneticRadius in subsequent calculations.

Lines 538--585: doing ODE integration by hand like this is both computationally costly and error-prone. I highly recommend using an off-the-shelf ODE integrator instead of using your own. See BaseBinaryStar::CalculateMassTransferOrbit() for an example.

There are some random 1000.0 hanging around; are these meant to represent g to kg conversion? if so, use G_TO_KG, will make the code clearer for a reader.

There is an int divideTimestepBy that's used to divide doubles; probably better to make this a double, to avoid giving the compiler freedom in (mis)interpreting type conversion in double/integer division -- but see above.

Not a fan of converting OPTIONS->NeutronStarAccretionInCE() into an integer NSCE and then using that integer; just use the option value directly, makes for cleaner code.

line 601:
double Jacc = MoI * PPOW(G_CGS * mass_g /r_cm_3, 0.5) * p_MassGainPerTimeStep * 1000.0 / mass_g ;
Why does the angular momentum of the accreted material care about the moment of inertia of the neutron star if you are assuming Keplerian accretion at the surface?

There's no need to define new variables like MoI, angularMomentum that just store values of existing variables -- only makes the code longer and harder to read.

[yuzhesong]

@ilyamandel working on cleaning up the code.

For your question regarding line 601:
double Jacc = MoI * PPOW(G_CGS * mass_g /r_cm_3, 0.5) * p_MassGainPerTimeStep * 1000.0 / mass_g ;

We are calculating the added angular momentum here as

J_acc = dJ = (J / M) * dm = MoI * sqrt(GM/R^3) * dm/M

====================
Also regarding your suggestion of using an off-the-shelf ODE integrator: I'm not particularly familiar with this tool, so I will have to look into implementing one for this part of the code. Do you have any other resources regarding this topic?

[ilyamandel]

Hi @yuzhesong :

General info on numerical ODE integration that uses higher-order integrators to improve accuracy/cost ratios and automatically chooses step sizes to achieve desired tolerances: Numerical Recipes textbook.

Specifically, I used Boost ODE integrators; there's documentation online (see https://www.boost.org/doc/libs/1_82_0/libs/numeric/odeint/doc/html/index.html ), but you are probably best off just using my example in BaseBinaryStar::CalculateMassTransferOrbit().

J_acc = dJ = (J / M) * dm = MoI * sqrt(GM/R^3) * dm/M

Are you assuming that the star is rotating at the Keplerian orbital frequency prior to accretion (that's what J = MoI * sqrt(GM/R^3) implies)?

[ilyamandel]

Closing at @yuzhesong 's suggestion since this is replaced by #1315 .