Update pulsar accretion treatment in NS.cpp

online-docs/pages/User guide/Program options/program-options-list-defaults.rst

@@ -845,6 +845,11 @@ Default = FIXED_MASS
 Amount of mass lost in neutrinos during BH formation (either as fraction or in solar masses, depending on the value of ``--neutrino-mass-loss-bh-formation``). |br|
 Default = 0.1
 
+**--neutron-star-accretion-in-ce** |br|
+Neutron star accretion behavior in common envelope. |br|
+Options: { ZERO, DISK, SURFACE } |br|
+Default = ZERO
+
 **--neutron-star-equation-of-state** |br|
 Neutron star equation of state. |br|
 Options: { SSE, ARP3 } |br|

src/BaseBinaryStar.cpp

@@ -388,11 +388,11 @@ void BaseBinaryStar::SetRemainingValues() {
     m_Flags.mergesInHubbleTime                       = false;
     m_Unbound                                        = false;
 
-    m_SystemicVelocity                           = Vector3d();
-    m_NormalizedOrbitalAngularMomentumVector     = Vector3d();
-	m_ThetaE                                     = DEFAULT_INITIAL_DOUBLE_VALUE;
-	m_PhiE                                       = DEFAULT_INITIAL_DOUBLE_VALUE;
-	m_PsiE                                       = DEFAULT_INITIAL_DOUBLE_VALUE;
+    m_SystemicVelocity                               = Vector3d();
+    m_NormalizedOrbitalAngularMomentumVector         = Vector3d();
+	m_ThetaE                                         = DEFAULT_INITIAL_DOUBLE_VALUE;
+	m_PhiE                                           = DEFAULT_INITIAL_DOUBLE_VALUE;
+	m_PsiE                                           = DEFAULT_INITIAL_DOUBLE_VALUE;
 
 	m_SynchronizationTimescale                       = DEFAULT_INITIAL_DOUBLE_VALUE;
 	m_CircularizationTimescale                       = DEFAULT_INITIAL_DOUBLE_VALUE;
@@ -1363,7 +1363,7 @@ bool BaseBinaryStar::ResolveSupernova() {
         Vector3d relativeVelocityVector       = relativeVelocityVectorPrev + (natalKickVector - companionRecoilVector);         // km/s - PostSN relative velocity vector
 
         Vector3d orbitalAngularMomentumVector = cross(separationVectorPrev, relativeVelocityVector);                            // km^2 s^-1 - PostSN specific orbital angular momentum vector
-        double   orbitalAngularMomentum = orbitalAngularMomentumVector.mag;                                                     // km^2 s^-1 - PostSN specific orbital angular momentum 
+        double   orbitalAngularMomentum       = orbitalAngularMomentumVector.mag;                                                     // km^2 s^-1 - PostSN specific orbital angular momentum 
         m_NormalizedOrbitalAngularMomentumVector = orbitalAngularMomentumVector/orbitalAngularMomentum;                         // set unit vector here to make printing out the inclination vector easier
 
         Vector3d eccentricityVector           = cross(relativeVelocityVector, orbitalAngularMomentumVector) / 
@@ -2166,8 +2166,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         }
     }
 
-	// Check for recycled pulsars. Not considering CEE as a way of recycling NSs.
-	if (!m_CEDetails.CEEnow && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) {                                           // accretor is a neutron star
+	// Check for recycled pulsars.  
+    if ((!m_CEDetails.CEEnow || OPTIONS->NeutronStarAccretionInCE() != NS_ACCRETION_IN_CE::ZERO)  && m_Accretor->IsOneOf({ STELLAR_TYPE::NEUTRON_STAR })) {
         m_Donor->SetRLOFOntoNS();                                                                                               // donor donated mass to a neutron star
         m_Accretor->SetRecycledNS();                                                                                            // accretor is (was) a recycled NS
 	}

src/NS.h

@@ -8,7 +8,7 @@
 
 #include "Remnants.h"
 #include "BH.h"
-
+#include <boost/numeric/odeint.hpp>
 
 class BaseStar;
 class Remnants;
@@ -46,8 +46,17 @@ class NS: virtual public BaseStar, public Remnants {
     static  double          CalculateRadiusOnPhase_Static(const double p_Mass)      { return CalculateRadiusOnPhaseInKM_Static(p_Mass) * KM_TO_RSOL; }              // Radius on phase in Rsol
 
     static  double          CalculateRemnantMass_Static(const double p_COCoreMass)  { return 1.17 + (0.09 * p_COCoreMass); }                                        // Hurley et al., eq 92
-
-
+    static  DBL_DBL_DBL_DBL  deltaAngularMomentumByPulsarAccretion(const double p_MassGainPerTimeStep, 
+                                                const double p_MagField, 
+                                                const double p_Mass, 
+                                                const double p_Radius,
+                                                const double p_SpinFrequency, 
+                                                const double p_AngularMomentum, 
+                                                const double p_Stepsize, 
+                                                const double p_Kappa, 
+                                                const double p_Epsilon,
+                                                const double p_MoI
+                                                ) ;
 protected:
 
     void Initialise() {
@@ -86,12 +95,22 @@ class NS: virtual public BaseStar, public Remnants {
 
     double          CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const;
 
-    double          ChooseTimestep(const double p_Time) const;
+    double          ChooseTimestep(const double p_Time) const;                              
+
+
 
     STELLAR_TYPE    EvolveToNextPhase()                                                     { return STELLAR_TYPE::BLACK_HOLE; }
 
     bool            ShouldEvolveOnPhase() const                                             { return (m_Mass <= OPTIONS->MaximumNeutronStarMass()); }               // Evolve as a neutron star unless mass > maximum neutron star mass (e.g. through accretion)
     void            SpinDownIsolatedPulsar(const double p_Stepsize);
+    DBL_DBL_DBL_DBL PulsarAccretion(const double p_MagField, 
+                                    const double p_SpinPeriod, 
+                                    const double p_AngularMomentum,
+                                    const double p_Stepsize,
+                                    const double p_MassGainPerTimeStep,
+                                    const double kappa,
+                                    const double p_Epsilon);
+
     void            UpdateMagneticFieldAndSpin(const bool   p_CommonEnvelope,
                                                const bool   p_RecycledNS,
                                                const double p_Stepsize,
@@ -100,4 +119,4 @@ class NS: virtual public BaseStar, public Remnants {
 
 };
 
-#endif // __NS_h__
+#endif // __NS_h__

src/Options.cpp

@@ -526,6 +526,9 @@ void Options::OptionValues::Initialise() {
     m_NeutronStarEquationOfState.type                               = NS_EOS::SSE;
     m_NeutronStarEquationOfState.typeString                         = NS_EOSLabel.at(m_NeutronStarEquationOfState.type);
 
+    // Neutron star accretion in common envelope
+    m_NSAccretionInCE.type                                          = NS_ACCRETION_IN_CE::ZERO;
+    m_NSAccretionInCE.typeString                                    = NS_ACCRETION_IN_CELabel.at(m_NSAccretionInCE.type) ; 
 
     // Pulsar birth magnetic field distribution
     m_PulsarBirthMagneticFieldDistribution.type                     = PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION::ZERO;
@@ -1807,7 +1810,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
             po::value<std::string>(&p_Options->m_NeutronStarEquationOfState.typeString)->default_value(p_Options->m_NeutronStarEquationOfState.typeString),                                                      
             ("Neutron star equation of state to use (" + AllowedOptionValuesFormatted("neutron-star-equation-of-state") + ", default = '" + p_Options->m_NeutronStarEquationOfState.typeString + "')").c_str()
         )
-
+        (
+            "neutron-star-accretion-in-CE",                              
+            po::value<std::string>(&p_Options->m_NSAccretionInCE.typeString)->default_value(p_Options->m_NSAccretionInCE.typeString),                                                      
+            ("Neutron star accretion in CE to use (" + AllowedOptionValuesFormatted("neutron-star-accretion-in-CE") + ", default = '" + p_Options->m_NSAccretionInCE.typeString + "')").c_str()
+        )
         (
             "OB-mass-loss",                                      
             po::value<std::string>(&p_Options->m_OBMassLoss.typeString)->default_value(p_Options->m_OBMassLoss.typeString),                                                                  
@@ -2247,6 +2254,11 @@ std::string Options::OptionValues::CheckAndSetOptions() {
             COMPLAIN_IF(!found, "Unknown Neutron Star Equation of State");
         }
 
+        if (!DEFAULTED("neutron-star-accretion-in-ce")) {                                                                         // neutron star accretion in ce
+            std::tie(found, m_NSAccretionInCE.type) = utils::GetMapKey(m_NSAccretionInCE.typeString, NS_ACCRETION_IN_CELabel, m_NSAccretionInCE.type);
+            COMPLAIN_IF(!found, "Unknown Neutron Star Accretion in CE");
+        }
+
         if (!DEFAULTED("OB-mass-loss")) {                                                                    // OB (main sequence) loss prescription
             std::tie(found, m_OBMassLoss.type) = utils::GetMapKey(m_OBMassLoss.typeString, OB_MASS_LOSS_LABEL, m_OBMassLoss.type);
             COMPLAIN_IF(!found, "Unknown OB Mass Loss Prescription");
@@ -2546,6 +2558,7 @@ std::vector<std::string> Options::AllowedOptionValues(const std::string p_Option
         case _("mode")                                              : POPULATE_RET(EVOLUTION_MODE_LABEL);                           break;
         case _("neutrino-mass-loss-BH-formation")                   : POPULATE_RET(NEUTRINO_MASS_LOSS_PRESCRIPTION_LABEL);          break;
         case _("neutron-star-equation-of-state")                    : POPULATE_RET(NS_EOSLabel);                                    break;
+        case _("neutron-star-accretion-in-CE")                      : POPULATE_RET(NS_ACCRETION_IN_CELabel);                                    break;
         case _("OB-mass-loss")                                      : POPULATE_RET(OB_MASS_LOSS_LABEL);                             break;
         case _("orbital-period-distribution")                       : POPULATE_RET(ORBITAL_PERIOD_DISTRIBUTION_LABEL);              break;
         case _("pulsar-birth-magnetic-field-distribution")          : POPULATE_RET(PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION_LABEL); break;
@@ -4624,6 +4637,7 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const {
         case PROGRAM_OPTION::NOTES                                          : value = Notes();                                                              break;
 
         case PROGRAM_OPTION::NS_EOS                                         : value = static_cast<int>(NeutronStarEquationOfState());                       break;
+        case PROGRAM_OPTION::NS_ACCRETION_IN_CE                             : value = static_cast<int>(NeutronStarAccretionInCE());                       break;
 
         case PROGRAM_OPTION::ORBITAL_PERIOD                                 : value = OrbitalPeriod();                                                      break;
         case PROGRAM_OPTION::ORBITAL_PERIOD_DISTRIBUTION                    : value = static_cast<int>(OrbitalPeriodDistribution());                        break;

src/Options.h

@@ -380,6 +380,8 @@ class Options {
         "maximum-mass-donor-nandez-ivanova",
         "minimum-secondary-mass",
 
+        "neutron-star-accretion-in-ce",
+
         "orbital-period",
         "orbital-period-distribution",
         "orbital-period-max",
@@ -510,6 +512,7 @@ class Options {
         "notes-hdrs",
         "neutrino-mass-loss-BH-formation",
         "neutron-star-equation-of-state",
+        "neutron-star-accretion-in-CE",
 
         "OB-mass-loss",
         "orbital-period-distribution",
@@ -965,6 +968,8 @@ class Options {
             // Neutron star equation of state
             ENUM_OPT<NS_EOS>                                    m_NeutronStarEquationOfState;                                   // NS EOS
 
+            // Neutron star accretion in common envelope
+            ENUM_OPT<NS_ACCRETION_IN_CE>                        m_NSAccretionInCE;                                             // NS Accretion In CE
 
             // Pulsar birth magnetic field distribution string
             ENUM_OPT<PULSAR_BIRTH_MAGNETIC_FIELD_DISTRIBUTION>  m_PulsarBirthMagneticFieldDistribution;                         // Birth magnetic field distribution for pulsars
@@ -1425,6 +1430,7 @@ class Options {
     double                                      NeutrinoMassLossValueBH() const                                         { return OPT_VALUE("neutrino-mass-loss-BH-formation-value", m_NeutrinoMassLossValueBH, true); }
 
     NS_EOS                                      NeutronStarEquationOfState() const                                      { return OPT_VALUE("neutron-star-equation-of-state", m_NeutronStarEquationOfState.type, true); }
+    NS_ACCRETION_IN_CE                          NeutronStarAccretionInCE() const                                        { return OPT_VALUE("neutron-star-accretion-in-ce", m_NSAccretionInCE.type, true); }
 
     std::string                                 Notes(const size_t p_Idx) const                                         { return OPT_VALUE("notes", m_Notes[p_Idx], true); }
     std::vector<std::string>                    Notes() const                                                           { return OPT_VALUE("notes", m_Notes, true); }

src/changelog.h

@@ -1225,8 +1225,16 @@
 //                                      - Change TPAGB::IsSupernova() so that stars with base of AGB core masses below MCBUR1 remain on the TPAGB until they make WDs; remove ResolveTypeIIaSN() functionality.
 //                                      - Add --evolve-main-sequence-mergers option which allows for main sequence merger products to continue evolution
 //                                      - Update HG::CalculateRadialExtentConvectiveEnvelope() to use a combination of Hurley & Picker to avoid excessively high convective envelope densities
+// 02.51.00    YS - Jul 03, 2024      - Update to neutron star accretion treatments:
+//                                      - Fixes to 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 programming option "NS-ACCRETION-IN-CE" for different treatment of how neutron star would behave when in CE. 
+
+
+const std::string VERSION_STRING = "02.51.00";
+
 
-const std::string VERSION_STRING = "02.50.00";
 
 
 # endif // __changelog_h__

src/compasConfigDefault.yaml

@@ -260,6 +260,7 @@ stringChoices:
 #    --fryer-supernova-engine: 'DELAYED'                                   # Default: 'DELAYED'                     # Options: ['DELAYED','RAPID']
 #    --kick-magnitude-distribution: 'MULLERMANDEL'                         # Default: 'MULLERMANDEL'                # Options: ['MULLERMANDEL','MULLER2016MAXWELLIAN','MULLER2016','BRAYELDRIDGE','MAXWELLIAN','FLAT','FIXED','ZERO']
 #    --kick-direction: 'ISOTROPIC'                                         # Default: 'ISOTROPIC'                   # Options: ['POLES','WEDGE','POWERLAW','PERPENDICULAR','INPLANE','ISOTROPIC']
+#    --neutron-star-accretion-in-ce: 'ZERO'                                # Default: 'ZERO'                        # Options: ['ZERO', 'DISK', 'SURFACE']
 #    --neutron-star-equation-of-state: 'SSE'                               # Default: 'SSE'                         # Options: ['ARP3','SSE']
 #    --neutrino-mass-loss-BH-formation: 'FIXED_MASS'                       # Default: 'FIXED_MASS'                  # Options: ['FIXED_MASS','FIXED_FRACTION']
 #    --pulsar-birth-magnetic-field-distribution: 'ZERO'                    # Default: 'ZERO'                        # Options: ['LOGNORMAL','UNIFORM','FLATINLOG','FIXED','ZERO']

src/constants.h

@@ -685,6 +685,7 @@ enum class ERROR: int {
     UNKNOWN_MT_THERMALLY_LIMITED_VARIATION,                         // unknown mass transfer thermally limited variation
     UNKNOWN_NEUTRINO_MASS_LOSS_PRESCRIPTION,                        // unknown neutrino mass loss prescription
     UNKNOWN_NS_EOS,                                                 // unknown NS equation-of-state
+    UNKNOWN_NS_ACCRETION_IN_CE,                                     // unknown NS accretion-in-ce
     UNKNOWN_PPI_PRESCRIPTION,                                       // unknown pulsational pair instability prescription
     UNKNOWN_PROGRAM_OPTION,                                         // unknown program option
     UNKNOWN_PROPERTY_TYPE,                                          // unknown property type
@@ -841,6 +842,7 @@ const COMPASUnorderedMap<ERROR, std::tuple<ERROR_SCOPE, std::string>> ERROR_CATA
     { ERROR::UNKNOWN_MASS_LOSS_PRESCRIPTION,                        { ERROR_SCOPE::ALWAYS,              "Unknown mass loss prescription" }},
     { ERROR::UNKNOWN_NEUTRINO_MASS_LOSS_PRESCRIPTION,               { ERROR_SCOPE::ALWAYS,              "Unknown neutrino mass loss prescription" }},
     { ERROR::UNKNOWN_NS_EOS,                                        { ERROR_SCOPE::ALWAYS,              "Unknown NS equation-of-state" }},
+    { ERROR::UNKNOWN_NS_ACCRETION_IN_CE,                            { ERROR_SCOPE::ALWAYS,              "Unknown NS accretion-in-ce" }},
     { ERROR::UNKNOWN_PPI_PRESCRIPTION,                              { ERROR_SCOPE::ALWAYS,              "Unknown pulsational pair instability prescription" }},
     { ERROR::UNKNOWN_PROGRAM_OPTION,                                { ERROR_SCOPE::ALWAYS,              "Unknown program option - property details not found" }},
     { ERROR::UNKNOWN_PROPERTY_TYPE,                                 { ERROR_SCOPE::ALWAYS,              "Unknown property type - property details not found" }},
@@ -1282,6 +1284,14 @@ const COMPASUnorderedMap<NS_EOS, std::string> NS_EOSLabel = {
     { NS_EOS::ARP3, "ARP3" }
 };
 
+// Neutron Star accretion in CE
+enum class NS_ACCRETION_IN_CE: int { ZERO, DISK, SURFACE };
+const COMPASUnorderedMap<NS_ACCRETION_IN_CE, std::string> NS_ACCRETION_IN_CELabel = {
+    { NS_ACCRETION_IN_CE::ZERO,     "ZERO" },
+    { NS_ACCRETION_IN_CE::DISK,     "DISK" },
+    { NS_ACCRETION_IN_CE::SURFACE,  "SURFACE" }
+};
+
 
 // Orbital Period Distributions
 enum class ORBITAL_PERIOD_DISTRIBUTION: int { FLATINLOG };
@@ -2587,6 +2597,7 @@ enum class PROGRAM_OPTION: int {
     NOTES,
 
     NS_EOS,
+    NS_ACCRETION_IN_CE,
 
     ORBITAL_PERIOD,
     ORBITAL_PERIOD_DISTRIBUTION,

src/typedefs.h

@@ -312,5 +312,6 @@ typedef struct StellarCEDetails {                      // Common Envelope detail
 typedef std::vector< double > state_type;
 typedef boost::numeric::odeint::runge_kutta_cash_karp54< state_type > error_stepper_type;
 typedef boost::numeric::odeint::controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
+typedef boost::numeric::odeint::runge_kutta4< state_type > stepper;
 
 #endif // __typedefs_h__

src/yaml.h

@@ -329,6 +329,7 @@ namespace yaml {
             "    --kick-magnitude-distribution",
             "    --kick-direction",
             "    --neutron-star-equation-of-state",
+            "    --neutron-star-accretion-in-ce",
             "    --neutrino-mass-loss-BH-formation",
             "    --pulsar-birth-magnetic-field-distribution",
             "    --pulsar-birth-spin-period-distribution",