Updated MSP treatments.

[yuzhesong]

Replaces #1082

Includes pulsar accretion treatments calculated with BOOST ODE integrator and other fixes, including solution to #1002

[ilyamandel]

Options.cpp:

    if (!DEFAULTED("neutron-star-accretion-in-ce")) {                                                                         // neutron star accretion in common envelope
        std::tie(found, m_NeutronStarAccretionInCE.type) = utils::GetMapKey(m_NeutronStarAccretionInCE.typeString, NS_ACCRETION_IN_CELabel, m_NeutronStarAccretionInCE.type);
        COMPLAIN_IF(!found, "Unknown Neutron Star Equation of State");
    }

Why is the complaint about an unknown NS EoS? This is about accretion, not EoS...

Same file,

    case _("neutron-star-accretion-in-ce")                      : POPULATE_RET(NS_ACCRETION_IN_CELabel);                        break;
    case _("neutron-star-equation-of-state")                    : POPULATE_RET(NS_EOSLabel);                                    break;

Format / capitalisation of labels different from other labels?

changelog.h:

• comments should be edited for English grammar
• "deltaAngularMomentumByPulsarAccretion()" -- that's not the name of the function, check capitalisation
• "- Update to NS::UpdateMagneticFieldAndSpin(), specifically adding" <- adding what?

A new option has been added. That means:

• the program options documentation needs to be updated accordingly
• a new default YAML should be created
• whats-new.rst should probably be updated

NS.cpp:

• comments for DeltaAngularMomentumByPulsarAccretion need to be adjusted for English grammar (e.g., "Calculate the magnetic field... and also returns"), correct capitalisation of function name, correct number and list of return variables in the tuple (the comment sometimes lists three, sometimes four -- the function returns three), correct ordering and list of input parameters, etc.

Do we need to pass parameters like p_Mass? This does not appear to be a static function, and a neutron star known its own mass -- would the argument passed in ever be different than m_Mass? Ditto for some of the other variables.

Why p_MassGainPerTimeStep? Isn't the important thing just the mass gain?

double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G;

If initialMagField_G was, say, zero to begin with, it will end up having a non-zero value but still one less than magFieldLowerLimit_G. I don't understand the logic...

More to come later.

[ilyamandel]

Aside: @yuzhesong , if this replaces #1082 , can we delete that one?

[yuzhesong]

@ilyamandel

Thanks for the review! Making changes to the codes...

To some of your points:

Do we need to pass parameters like p_Mass? This does not appear to be a static function, and a neutron star known its own mass -- would the argument passed in ever be different than m_Mass? Ditto for some of the other variables.

NS::DeltaAngularMomentumByPulsarAccretion() is used for the integration for NS in RLOF as a dynamic ODE, so certain parameters needs to be updated. mass, magnetic field, spin frequency etc are all updated throughout all iterations.

Why p_MassGainPerTimeStep? Isn't the important thing just the mass gain?

p_MassGainPerTimeStep is the mass gain for that specific time step. If the name is confusing, I can change it to p_MassGain.

double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G;

If initialMagField_G was, say, zero to begin with, it will end up having a non-zero value but still one less than magFieldLowerLimit_G. I don't understand the logic...

I don't quite understand the issue with this? It's following Eq. 12 in arxiv:1912.02415 to calculate the decay of magnetic field from mass accretion. There wouldn't be any physical situation in a NS's evolution for the magnetic field to be low enough, certainly not zero. The magnetic field needs to be evolved from a certain value, either decays through spin down or mass accretion. The birth value is drawn from a distribution specified by the program options, and the lower limit of magnetic field should be always lower than that entire distribution.

Aside: @yuzhesong , if this replaces #1082 , can we delete that one?

Yes, #1082 can be deleted

[jeffriley]

Placeholder - I have some changes I'd like made.

[ilyamandel]

Just a quick note that I have to go offline now (else I'll have to swim over to Tasmania as the ferry will leave without me), so I'll have to return to the review, which I haven't quite finished, in January of 2025. Happy New Year!

[jeffriley]

In addition to @ilyamandel's comments:

In void NS::CalculateAndSetPulsarParameters() you have the following code:

m_PulsarDetails.spinPeriod    = CalculateBirthSpinPeriod();                                         // spin period in ms
m_PulsarDetails.spinFrequency = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);

CalculateBirthSpinPeriod() can return a value of 0.0 by design (e.g. when OPTIONS->PulsarBirthSpinPeriodDistribution() == ZERO, which is actually the default), which results in a divide-by-zero error - and sets
m_PulsarDetails.spinFrequency to INF, and we then use that value as a parameter to NS::CalculateSpinDownRate(), propagating the error...

At a minimum there should be a check here for spin period = 0.0, but if the user has specified pulsar-birth-spin-period-distribution ZERO (or left it at the default value), it seems to me we don't need to call NS::CalculateAndSetPulsarParameters() at all,

NS::SpinDownIsolatedPulsar() also divides by m_PulsarDetails.spinFrequency - do we even need to call NS::SpinDownIsolatedPulsar() if m_PulsarDetails.spinFrequency == 0.0? At least put a check at the top of the function to check for zero spin and exit immediately.

Ditto for NS::UpdateMagneticFieldAndSpin() - can wet put a check at the top of the function to check for zero spin and exit immediately?

How much of the pulsar code can we actually avoid (and should avoid to prevent divide-by-zero errors) if we know the star is not spinning?

[yuzhesong]

@jeffriley

Added solution to issue #1257

[ilyamandel]

@yuzhesong :

p_MassGainPerTimeStep is the mass gain for that specific time step. If the name is confusing, I can change it to p_MassGain.

That sounds like a good solution.

double newPulsarMagneticField = (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G;

If initialMagField_G was, say, zero to begin with, it will end up having a non-zero value but still one less than magFieldLowerLimit_G. I don't understand the logic...

I don't quite understand the issue with this? It's following Eq. 12 in arxiv:1912.02415 to calculate the decay of magnetic field from mass accretion. There wouldn't be any physical situation in a NS's evolution for the magnetic field to be low enough, certainly not zero. The magnetic field needs to be evolved from a certain value, either decays through spin down or mass accretion. The birth value is drawn from a distribution specified by the program options, and the lower limit of magnetic field should be always lower than that entire distribution.

I would not assume that initialMagField_G is lower than magFieldLowerLimit_G. You never know in advance what silly choices someone else might make. :-) A safer solution could be as follows:

double newPulsarMagneticField = initialMagField_G < magFieldLowerLimit_G ? magFieldLowerLimit_G : (initialMagField_G - magFieldLowerLimit_G) * exp(-p_MassGainPerTimeStep / p_Kappa) + magFieldLowerLimit_G;

[ilyamandel]

double p = R_CM_6 * R_CM_6 / (mass_g * mass_g * mDot);

in NS::DeltaAngularMomentumByPulsarAccretion_Static()

This is not correct and does not match Eq. (10) of https://www.arxiv.org/pdf/1912.02415 (see powers of M and Mdot).

@yuzhesong , typos happen, but this is going very slowly because typos seem to be particularly common in this code and we are all busy... I suggest that you check that your code is yielding correct results (as compared to some alternative code that you write independently, say, in Python or Matlab) for a range of parameter choices -- otherwise, we'll be reviewing this for a while. ;-)

[Aside: if you haven't done it yet, I would also personally recommend re-deriving quantities like the Alfven radius yourself; I find that once you've derived something and therefore understood the scalings, copy-paste typo errors when copying equations are less likely.]

[ilyamandel]

Is it not possible the value on entry to the function is set?

If the frequency is set, we should not be calling
m_PulsarDetails.spinPeriod = CalculateBirthSpinPeriod();
-- we should be using the frequency to set the period. The current behaviour does not appear to be self-consistent or to match the desired goals.

[ilyamandel]

EDIT: I did test this when I implemented the error handling code (that works on Linux...), and dividing a floating-point value by 0.0 sets +/- inf appropriately. I can't remember if I tested dividing a floating-point valued by integer 0.

But why not just use the following (we already do in various other places in the code) instead of trying to divide by zero? Seems safer...
std::numeric_limits::infinity();

[yuzhesong]

@ilyamandel

I also still don't understand (and didn't see a response above to my previous change request on this) why a magnetic field of zero forces a change in the spin period.

Essentially we consider if a neutron star with zero spin frequency, infinite spin period, or zero magnetic field to be non-spinning in this situation, so we force all the spin-related parameters to return 0 for the neutron star under this situation. This is accounting for the option of setting magnetic field at birth to ZERO in the program option.

[jeffriley]

I agree, which is why I suggested it in other places. I think Robert wanted to emphasise the relationship, but upon reflection I think setting it explicitly is better.

…

On Tue, 11 Mar 2025, 9:28 pm Ilya Mandel, ***@***.***> wrote: EDIT: I did test this when I implemented the error handling code (that works on Linux...), and dividing a floating-point value by 0.0 sets +/- inf appropriately. I can't remember if I tested dividing a floating-point valued by integer 0. But why not just use the following (we already do in various other places in the code) instead of trying to divide by zero? Seems safer... std::numeric_limits::infinity(); — Reply to this email directly, view it on GitHub <#1315 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGM5XNROUUWB5J7TYR3SMSD2T23DLAVCNFSM6AAAAABT3ZWWK6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMJTGU2TSNZUGI> . You are receiving this because you were mentioned.Message ID: ***@***.***> [image: ilyamandel]*ilyamandel* left a comment (TeamCOMPAS/COMPAS#1315) <#1315 (comment)> EDIT: I did test this when I implemented the error handling code (that works on Linux...), and dividing a floating-point value by 0.0 sets +/- inf appropriately. I can't remember if I tested dividing a floating-point valued by integer 0. But why not just use the following (we already do in various other places in the code) instead of trying to divide by zero? Seems safer... std::numeric_limits::infinity(); — Reply to this email directly, view it on GitHub <#1315 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGM5XNROUUWB5J7TYR3SMSD2T23DLAVCNFSM6AAAAABT3ZWWK6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMJTGU2TSNZUGI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[yuzhesong]

EDIT: I did test this when I implemented the error handling code (that works on Linux...), and dividing a floating-point value by 0.0 sets +/- inf appropriately. I can't remember if I tested dividing a floating-point valued by integer 0.

But why not just use the following (we already do in various other places in the code) instead of trying to divide by zero? Seems safer... std::numeric_limits::infinity();

@ilyamandel i'm happy with doing it that way.

[ilyamandel]

NS.cpp, lines 448, 449:

double omegaK = std::sqrt(G_CGS * p_Mass / magneticRadius) / magneticRadius;                // Hz
double vDiff  = omegaK - p_SpinFrequency;                                                   // Hz

Although Hz and rad/s are dimensionally equivalent, by convention, angular velocities are always quoted in rad/s, Hz is reserved for frequencies.

===========

Line 508:

    double jAcc            = m_MomentOfInertia_CGS * std::sqrt(G_CGS * massG / (radiusCM * radiusCM * radiusCM)) * p_MassGain / G_TO_KG / massG;    
    m_AngularMomentum_CGS += jAcc;  

I didn't quite understand what SURFACE means, but it looks like it just means "keep angular velocity fixed". In that case, what's the point of computing fDot on line 520 -- should it not be identically zero?

=====

[ilyamandel]

@ilyamandel

I also still don't understand (and didn't see a response above to my previous change request on this) why a magnetic field of zero forces a change in the spin period.

Essentially we consider if a neutron star with zero spin frequency, infinite spin period, or zero magnetic field to be non-spinning in this situation, so we force all the spin-related parameters to return 0 for the neutron star under this situation. This is accounting for the option of setting magnetic field at birth to ZERO in the program option.

That's not accurate: ZERO and 0.0 are not the same here. I already commented on this in a previous review.

    case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO:                                                    // ZERO
        log10B = 0.0;
        break;

This sets the log of the B field to zero. That means the B field is set to be 1.0 Gauss, not 0.0.

Meanwhile, the code in question looks at the comparison

utils::Compare(m_PulsarDetails.magneticField, 0.0) .

1.0 \neq 0.0, so with the option ZERO, this will in fact not match.

[yuzhesong]

@ilyamandel

I also still don't understand (and didn't see a response above to my previous change request on this) why a magnetic field of zero forces a change in the spin period.

Essentially we consider if a neutron star with zero spin frequency, infinite spin period, or zero magnetic field to be non-spinning in this situation, so we force all the spin-related parameters to return 0 for the neutron star under this situation. This is accounting for the option of setting magnetic field at birth to ZERO in the program option.

That's not accurate: ZERO and 0.0 are not the same here. I already commented on this in a previous review.

    case PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO:                                                    // ZERO
        log10B = 0.0;
        break;

This sets the log of the B field to zero. That means the B field is set to be 1.0 Gauss, not 0.0.

Meanwhile, the code in question looks at the comparison

utils::Compare(m_PulsarDetails.magneticField, 0.0) .

1.0 \neq 0.0, so with the option ZERO, this will in fact not match.

Ah you're right. The magnetic field distribution setting to ZERO should return ZERO magnetic field not 1.0. I will change this.

[yuzhesong]

Just a note: @SimonStevenson and I were discussing this earlier today, and we thought that these unphysical behaviours (0 spin frequency, infinity spin period or 0 magnetic field) really shouldn't exist, because why record pulsar evolution into a file when the pulsars are not evolving...

Would you @ilyamandel and @jeffriley oppose to remove these setups? It would actually automatically nullify these discussions. We can still have the checks in the calculations just in case something goes wrong but we really should not encourage the setting of non-spinning pulsars here.

[jeffriley]

we really should not encourage the setting of non-spinning pulsars here.

I agree with that. If we provide the "NOSPIN" option, aren't we really just saying we will treat all neutron stars as vanilla neutron stars, and none as pulsars (effectively turning off pulsar evolution in COMPAS)? Isn't that the same (or should it be the same) as the --evolve-pulsars option? Also, zero spin frequency is unphysical for all stars isn't it (maybe not in COMPAS...)?

Here's my naive take on it. All pulsars are neutron stars, but not all neutron stars are pulsars - that's right, isn't it? Maybe we should have implemented pulsars as their own class, inheriting from the NS class - and maybe even their own stellar (sub?)type. If we actually have a pulsar class (and maybe stellar type), then we could define the spin threshold at which a NS spins up and becomes a pulsar (presumably until they spin fast enough there are no jets?), and the conditions under which it spins down and there are no longer jets and it is no longer a pulsar (should it then just become a NS again?). Should we consider a rewrite (really shouldn't be that hard), or is that something for the future?

[yuzhesong]

All pulsars are neutron stars, but not all neutron stars are pulsars - that's right, isn't it?

This is correct. But, I think in practice (or observation) the differentiation between a NS and a pulsar is if it's observed to pulsate, and there are many reasons why the pulsation is not observed (beaming fraction, beaming angle and direction, spin down luminosity is too low, too far away, to name a few). Plus, in theory, due to the conservation of angular momentum, all neutron stars born from supernovae should have a spin period and it should be relatively low (fast spinning at birth).

then we could define the spin threshold at which a NS spins up and becomes a pulsar (presumably until they spin fast enough there are no jets?), and the conditions under which it spins down and there are no longer jets and it is no longer a pulsar (should it then just become a NS again?).

So for all that above I personally would just recommend us getting rid of the no spin and zero magnetic field options, instead of creating an additional class just to accommodate the no-spinning cases. If people don't need the pulsar informations, they don't need to set --evolve-pulsar=True. All other stellar/binary properties of the neutron star systems (mass, radius, orbits, mass transfer, etc) should already be recorded even if --evolve-pulsar=False. But there are chances I might be missing something here? @SimonStevenson correct me here if I'm wrong.

[jeffriley]

instead of creating an additional class just to accommodate the no-spinning cases

No, maybe I didn't explain clearly. Pulsars are types of neutron stars aren't they. That lends itself to having a pulsar class that inherits from the NS class. Then all the code that relates specifically to pulsars, but not to plain neutron stars, resides in the pulsar class - and then, presumably, we can make assumptions in the pulsar class because we wouldn't be there if those assumptions were not true. And then when the pulsars spin down they go back to being plain neutron stars? And we don't have to worry about the checks in the pulsar code.

All other stellar/binary properties of the neutron star systems (mass, radius, orbits, mass transfer, etc) should already be
recorded even if --evolve-pulsar=False.

That's right - that would be the NS class, as opoosed to the Pulsar class.

If people don't need the pulsar informations, they don't need to set --evolve-pulsar=True.

True, and that was my point - or question. What's the difference between --evolve-pulsars False and --evolve-pulsars True + --pulsar-birth-spin-period-distribution NOSPIN - aren't they effectively the same thing? Or is that simplifying too much?

[yuzhesong]

I think I get your point now. But -

What's the difference between --evolve-pulsars False and --evolve-pulsars True + --pulsar-birth-spin-period-distribution NOSPIN - aren't they effectively the same thing? Or is that simplifying too much?

Logically yes they are the same thing. What the code is doing right now is that:

• --evolve-pulsars False means we end up with a neutron star and that's it. It doesn't call on anything in NS.cpp to calculate spin/magnetic field or anything.
• --pulsar-birth-spin-period-distribution NOSPIN creates the SSE/BSE_Pulsar_Evolution file and records everything (except for mass/radii/orbital properties) to be zero. Which really is unnecessary.

All other stellar/binary properties of the neutron star systems (mass, radius, orbits, mass transfer, etc) should already be
recorded even if --evolve-pulsar=False.

That's right - that would be the NS class, as opoosed to the Pulsar class.

These are in the BaseStar or BaseBinaryStar classes, so maybe not necessary to have an additional class just for them when the star is a NS?

So I guess my point is that we already have all the necessary parameters for a neutron star calculated and stored as stellar properties. Do we really need another class for NSs that are not pulsars?

[jeffriley]

Logically yes they are the same thing

So that's my point - why would anybody enable pulsars and the set the birth spin period to zero? They wouldn't, right? So we should remove NOSPIN.

These are in the BaseStar or BaseBinaryStar classes

All (or most) stellar properties are recorded in the base class, and are inherited by the derived classes.

, so maybe not necessary to have an additional class just for them when the star is a NS?

I'm not sure I understand that comment - we already have a class for when a star is a NS - it's the NS class...

Do we really need another class for NSs that are not pulsars?

No, that's not what I'm saying. We already have the NS class - that will stay, and that is for all neutron stars. The (new) pulsar class would be for neutron stars that are also pulsars. The Pulsar class would inherit from the NS class - it's just an OO way of looking at, and implementing, it...

[jeffriley]

Ok, so @yuzhesong has convinced me - either all neutron stars in COMPAS are pulsars, or none are - there is no physical distinction that would allow us to decide if/when a NS is actually a pulsar or not. So no need for a separate pulsar class.

[yuzhesong]

Line 508:

    double jAcc            = m_MomentOfInertia_CGS * std::sqrt(G_CGS * massG / (radiusCM * radiusCM * radiusCM)) * p_MassGain / G_TO_KG / massG;    
    m_AngularMomentum_CGS += jAcc;  

I didn't quite understand what SURFACE means, but it looks like it just means "keep angular velocity fixed". In that case, what's the point of computing fDot on line 520 -- should it not be identically zero?

@ilyamandel SURFACE here means pulsar accretes material directly onto its surface from the common envelope, so the jAcc here calculates the increase in angular momentum by adding p_MassGain mass at the co-rotation angular velocity at NS surface, while the NS is spinning at its own spin frequency. So fDot here shouldn't be zero.

Other than that, following the discussion above, I have deprecated the "ZERO" option for --pulsar-birth-magnetic-field-distribution and "NOSPIN" (originally "ZERO") option for --pulsar-birth-spin-period-distribution. I also removed all the guards in NS.cpp as there is no need for them if we're not dealing with non spinning NSs any more.

[ilyamandel]

@ilyamandel SURFACE here means pulsar accretes material directly onto its surface from the common envelope, so the jAcc here calculates the increase in angular momentum by adding p_MassGain mass at the co-rotation angular velocity at NS surface, while the NS is spinning at its own spin frequency. So fDot here shouldn't be zero.

Just to repeat what I said on the call: in that case, line 475 of NS.cpp is not correct. The specific Keplerian angular momentum at the surface is just sqrt(GM_NSR_NS), and j_Acc should be p_MassGain times this quantity. It should not depend on the NS moment of inertia.

Other than that, following the discussion above, I have deprecated the "ZERO" option for --pulsar-birth-magnetic-field-distribution and "NOSPIN" (originally "ZERO") option for --pulsar-birth-spin-period-distribution. I also removed all the guards in NS.cpp as there is no need for them if we're not dealing with non spinning NSs any more.

OK!

But you accept a draw of a zero spin period; if you really don't want to allow them (I know, set of measure zero, but still), you'll want to change < to <= on line 167 on NS.cpp:

do { pSpin = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_SPIN_ITERATIONS && utils::Compare(pSpin, 0.0) < 0);

to

do { pSpin = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_SPIN_ITERATIONS && utils::Compare(pSpin, 0.0) <= 0);

[ilyamandel]

Not sure we really need this around line 225 of NS.cpp: even if lof10B is negative, the actual B field is always positive...

        // this should terminate naturally, but just in case we add a guard
        std::size_t iterations = 0;
        do { log10B = RAND->RandomGaussian(sigma) + mean;} while (iterations++ < PULSAR_MAG_ITERATIONS && utils::Compare(log10B, 0.0) < 0);
        if (iterations >= PULSAR_MAG_ITERATIONS) THROW_ERROR(ERROR::TOO_MANY_PULSAR_MAG_ITERATIONS);

Would it be more useful to replace
utils::Compare(log10B, 0.0) < 0
with utils::Compare(log10B, log10(NS::NS_MAG_FIELD_LOWER_LIMIT) < 0
to ensure that the B field is always set to be above NS::NS_MAG_FIELD_LOWER_LIMIT ?

[ilyamandel]

The frequency derivative is not computed correctly on line 487 of NS.cpp:

double fDot = jAcc / p_Stepsize / m_MomentOfInertia_CGS;

There are two reasons why this isn't right. First, the change in frequency is not the change in angular momentum divided by the moment of inertia. One can add angular momentum and keep the frequency the same as it was. Secondly, this doesn't account for the change in moment of inertia (which may have changed because mass was added).

Fortunately, there's an easy way to do the calculation correctly. Replace the following lines

m_PulsarDetails.spinFrequency = m_AngularMomentum_CGS / m_MomentOfInertia_CGS;
m_PulsarDetails.spinPeriod = _2_PI / m_PulsarDetails.spinFrequency;
double fDot = jAcc / p_Stepsize / m_MomentOfInertia_CGS;

with

double previousSpinFrequency = m_PulsarDetails.spinFrequency;
m_PulsarDetails.spinFrequency = m_AngularMomentum_CGS / m_MomentOfInertia_CGS;
m_PulsarDetails.spinPeriod = _2_PI / m_PulsarDetails.spinFrequency;
double fDot = (m_PulsarDetails.spinFrequency - previousSpinFrequency) / p_Stepsize;

[ilyamandel]

Woohoo, we are done! ;-)