Magnetic field decay of NSs due to mass accretion

[astro-sjgao]

I am using COMPAS to calculate the pulsar evolution in X-ray binaries. However, I have observed that the magnetic field decay does not seem to be working as expected.

I have identified a potential cause in the file NS.cpp, as shown in the following line:
double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-1 * p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit ;

p_MassGainPerTimeStep seems to be in units of kg (as I found in BinaryConstituentStar.h), while kappa is also in units of kg. The term /1000.0 seems to be a typo.

So, the magnetic filed dacay resulting from mass accretion is insignificant.

For example, assuming the massscale is 0.02 solar mass, and a pulsar accreted 0.01 solar mass, the magnetic field B (>>Bmin) should be reduced to
(B-Bmin)*exp(-0.01/0.02)+Bmin~B*exp(-0.5)=0.61*B
rather than
(B-Bmin)*exp(-0.01/1000/0.02)+Bmin~B*exp(-0.5/1000)=0.9995*B.

[ilyamandel]

Dear ShiJie Gao, Thank you for your question. I am cc’ing Simon and Debatri, who should be able to help you out. Best wishes, Ilya

…

On 20 Oct 2023, at 09:47, ShiJie Gao ***@***.***> wrote: I am using COMPAS to calculate the pulsar evolution in X-ray binaries. However, I have observed that the magnetic field decay does not seem to be working as expected. I have identified a potential cause in the file NS.cpp, as shown in the following line: double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-1 * p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit ; p_MassGainPerTimeStep seems to be in units of kg (as I found in BinaryConstituentStar.h), while kappa is also in units of kg. So, the magnetic filed dacay resulting from mass accretion is insignificant. For example, assuming the massscale is 0.02 solar mass, and a pulsar accreted 0.01 solar mass, the magnetic field B (>>Bmin) should be reduced to (B-Bminexp(-0.01/0.02)+Bmin~Bexp(-0.5)=0.61B rather than (B-Bmin)exp(-0.01/1000/0.02)+Bmin~Bexp(-0.5/1000)=0.9995B. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[SimonStevenson]

Hi Shijie,

I believe that you have been in contact with @yuzhesong over email about this issue.

I haven't really looked at this version of NS.cpp, but it does seem to have a few issues. I think that you are right that the factor of 1000 must be a typo (even if it was meant to be a unit conversion, this wouldn't be the right way to do it). It is possible that someone thought that the equation needed to be in cgs units (so masses in grams) in which case there is a factor of 1000 in the conversion (but we should use something like KG_TO_G or MSOL_TO_G to convert). However, in this case, this factor would need to be applied to both the numerator and denominator of this fraction (both have to be in the same units). At present, it appears the numerator is in grams, and the denominator is in kg.

Typically in COMPAS, masses are in solar masses. So it is already a bit weird that p_MassGainPerTimeStep is in kg.

There are a couple of places where we have this rogue factor of 1000:

        double newPulsarMagneticField = (initialMagField - magFieldLowerLimit) * exp(-1 * p_MassGainPerTimeStep / 1000.0 / kappa) + magFieldLowerLimit ;

and later we have

       double mDot          =  p_MassGainPerTimeStep / 1000.0 / p_Stepsize ;

There seems to be a bit of a nasty mix of SI and cgs units throughout NS.cpp. It seems that similar magic factors appear elsewhere in NS.cpp: e.g., a magic factor of 1000000 appears in the line

    constexpr double _3_C_3         = 3.0 * C * C * C * 1000000.0;

in SpinDownIsolatedPulsar (I think this is the conversion from m to cm, cubed). I think we need to go through and clean this file up. I can try to do that, and put together a PR to address some of these changes, but I think a larger overhaul is needed here.

As @yuzhesong mentioned, although accretion should be possible through both mass transfer and common envelope evolution, in Chattopadhyay+2020 (https://arxiv.org/abs/1912.02415) we focused on accretion during the common envelope. Certainly for MSP formation, accretion through stable mass transfer will be more important. We have yet to explore this with COMPAS.

I also find it a bit weird that UpdateMagneticFieldAndSpin takes 5 parameters in NS.cpp but is called with 3 parameters in BaseBinaryStar.cpp. This seems confusing (at best), and possibly wrong (at worst).

TLDR: Yes, I believe the factor of 1000 is a left over/unneeded cgs conversion. This function should be considered experimental/unsupported at present.

[astro-sjgao]

Dear Ilya and Simon, thanks for your kind replies.

I have two additional comments on the plusar recycling process.

• super-Eddingtion accretion during CEE

  In chattopadhyay+2020, mass acccretion by an NS is based on the MacLeod's model, which assumes hyper-Eddington accretion to ensure that NS accretes enough material.

  It appears that COMPAS does not model the mass loss process during CEE step by step in detail, but instead directly removes the entire envelope.

  So, the assumed timestep during CEE determinates the average mass accretion rate of NSs, and the calculated mass accretion rate could be super-Eddington in some cases without a doubt. However, the mass accretion rate still seems to be limited to the Eddington accretion rate in Chattopadhyay+2020, (see Equation 16 in their paper), which contradicts the MacLeod's model.

• Formation of milisecond pulsars via RLOF

  During low-mass X-ray binary stage, the mass donor can be a main-sequence star without a clear or even without a helium core (Case A mass transfer). This recycling process may differ from the recycling process during CEE.

Best wishes,
Shijie, Gao

[astro-sjgao]

Dear Ilya, Simon and Yuzhe,

I found a typo in Equation 9 of Chattopadhyay+2020 paper, where the right-hand side is missing a factor of 2. The correct form of the equation should be as follows:

$$\frac{1}{\Omega_{\rm f}^2}-\frac{1}{\Omega_{\rm i}^2}=\frac{P_{\rm f}^2}{4\pi^2}-\frac{P_{\rm i}^2}{{4}\pi^2}=\left[\frac{4\pi}{\mu_0}\right]\frac{\textcolor{red}{4} R^6\sin^2 \alpha}{3c^3I}\left[B_{\rm min}^2(t_{\rm f}-t_{\rm i})-\tau B_{\rm min}(B_{\rm f}-B_{\rm i})-\frac{\tau}{2}(B_{\rm f}^2-B_{\rm i}^2)\right].$$

See also Equation 9 in Sgalletta+2023 (https://doi.org/10.1093/mnras/stad2768). There is another typo in their Equation 9, where the term "P" is unnecessary and can be removed.

Additionally, I checked the corresponding code in NS.cpp (lines 384-389)and found that there is also a missing factor of 2.

Best wishes,
Shijie, Gao

Formula derivation:

[debatric]

Hi Shijie, I am in a ship to Antarctica at the moment with very choppy internet. The typo in the paper was too late to correct, the correct equations were implemented in the code and my thesis which is available online. Regards, Debatri.

…

On Thu, Nov 16, 2023 at 8:54 AM ShiJie Gao ***@***.***> wrote: Dear Ilya, Simon and Yuzhe, I found a typo in Equation 9 of Chattopadhyay+2020 paper, where the right-hand side is missing a factor of 2. The correct form of the equation should be as follows: $$\frac{1}{\Omega_{\rm f}^2}-\frac{1}{\Omega_{\rm i}^2}=\frac{P_{\rm f}^2}{4\pi^2}-\frac{P_{\rm i}^2}{{4}\pi^2}=\left[\frac{4\pi}{\mu_0}\right]\frac{\textcolor{red}{4} R^6\sin^2 \alpha}{3c^3I}\left[B_{\rm min}^2(t_{\rm f}-t_{\rm i})-\tau B_{\rm min}(B_{\rm f}-B_{\rm i})-\frac{\tau}{2}(B_{\rm f}^2-B_{\rm i}^2)\right].$$ See also Equation 9 in Sgalletta+2023 ( https://doi.org/10.1093/mnras/stad2768). There is another typo in their Equation 9, where the term "P" is unnecessary and can be removed. Additionally, I checked the corresponding code in NS.cpp (lines 384-389)and found that there is also a missing factor of 2. Best wishes, Shijie, Gao Formula derivation: [image: 截屏2023-11-16 19 52 46] <https://user-images.githubusercontent.com/70212620/283439700-b24cc35c-e612-4720-8f63-a8571a63c8b3.png> — Reply to this email directly, view it on GitHub <#1002 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AI2EB52EIXUBZX27IUSEP7TYEX5IHAVCNFSM6AAAAAA6ITJU4KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMJUGI4TSMZXGE> . You are receiving this because you were assigned.Message ID: ***@***.***>

[ilyamandel]

@debatric , @SimonStevenson , @yuzhesong -- is this an issue on its own, separate from the issue described in #1082 ? [I note that Debatri wrote "The typo in the paper was too late to correct, the correct equations were implemented in the code and my thesis which is available online" -- if so, is it just the case of updating the comments in the code to point to the correct equations? I'd like to resolve long-standing issues, especially with bug labels...

[yuzhesong]

@ilyamandel #1082 contains the solutions to this issue. Just need to get it to the stage of being able to merged back into dev.