Implementing the Shikauchi+ (2024) core mass prescription for main sequence stars

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -236,6 +236,7 @@ stringChoices:
 #    --initial-mass-function: 'KROUPA'                                     # Default: 'KROUPA'                      # Options: ['KROUPA','UNIFORM','POWERLAW','SALPETER']
 #    --LBV-mass-loss-prescription: 'HURLEY_ADD'                            # Default: 'HURLEY_ADD'                  # Options: ['BELCZYNSKI','HURLEY','HURLEY_ADD','ZERO','NONE']
 #    --luminous-blue-variable-prescription: 'HURLEY_ADD'                   # Default: 'HURLEY_ADD'                  # Options: ['BELCZYNSKI','HURLEY','HURLEY_ADD','ZERO','NONE']
+#    --main-sequence-core-mass-prescription: 'MANDEL'                      # Default: 'MANDEL'                      # Options: ['SHIKAUCHI','MANDEL','NONE']
 #    --mass-loss-prescription: 'MERRITT2024'                               # Default: 'MERRITT2024'                 # Options: ['MERRITT2024','BELCZYNSKI2010','HURLEY','ZERO','NONE']
 #    --OB-mass-loss: 'VINK2021'                                            # Default: 'VINK2021'                    # Options: ['KRTICKA2018','BJORKLUND2022','VINK2021','VINK2001','ZERO','NONE']
 #    --OB-mass-loss-prescription: 'VINK2021'                               # Default: 'VINK2021'                    # Options: ['KRTICKA2018','BJORKLUND2022','VINK2021','VINK2001','ZERO','NONE']

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

@@ -795,6 +795,15 @@ Default = 4.2
 
 :ref:`Back to Top <options-props-top>`
 
+**--main-sequence-core-mass-prescription** |br|
+Main sequence core mass prescription. |br|
+Options: {NONE, MANDEL, SHIKAUCHI} |br|
+``NONE``      : No core mass treatment |br|
+``MANDEL``    : The core following case A mass transfer is set equal to the expected core mass of a newly formed HG star
+ with mass equal to that of the donor, scaled by the fraction of the donor's MS lifetime at mass transfer |br|
+``SHIKAUCHI`` : Core mass according to Shikauchi et al. (2024) |br|
+Default = MANDEL |br|
+
 **--mass-change-fraction** |br|
 Approximate desired fractional change in stellar mass on phase when setting SSE and BSE timesteps (applied before ``--timestep--multiplier``). |br|
 Recommended value is 0.005. |br|
@@ -1161,7 +1170,8 @@ Default = MULLERMANDEL
 If TRUE, preserve a larger donor core mass following case A mass transfer. |br|
 The core is set equal to the expected core mass of a newly formed HG star with mass equal to that of the donor,
 scaled by the fraction of the donor's MS lifetime at mass transfer. |br|
-Default = TRUE
+Default = TRUE |br|
+DEPRECATION NOTICE: this option has been deprecated and will soon be removed. Please use ``--main-sequence-core-mass-prescription MANDEL`` in future.
 
 **--revised-energy-formalism-nandez-ivanova** |br|
 Enable revised energy formalism of Nandez & Ivanova. |br|
@@ -1418,7 +1428,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 **Stellar evolution and winds**
 
 --use-mass-loss, --check-photon-tiring-limit, --cool-wind-mass-loss-multiplier, --luminous-blue-variable-prescription, --LBV-mass-loss-prescription
---luminous-blue-variable-multiplier, --mass-loss-prescription, --overall-wind-mass-loss-multiplier, --wolf-rayet-multiplier, 
+--luminous-blue-variable-multiplier, --main-sequence-core-mass-prescription, --mass-loss-prescription, --overall-wind-mass-loss-multiplier, --wolf-rayet-multiplier, 
 --expel-convective-envelope-above-luminosity-threshold, --luminosity-to-mass-threshold,
 --OB-mass-loss, --OB-mass-loss-prescription, --RSG-mass-loss, --RSG-mass-loss-prescription, --VMS-mass-loss, --vms-mass-loss-prescription, --WR-mass-loss, --WR-mass-loss-prescription
 

src/BaseBinaryStar.cpp

@@ -1730,6 +1730,16 @@ void BaseBinaryStar::ResolveMainSequenceMerger() {
 
     m_SemiMajorAxis = std::numeric_limits<float>::infinity();                                   // set separation to infinity to avoid subsequent fake interactions with a massless companion (RLOF, CE, etc.)
 
+    if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_Star1->MZAMS(), SHIKAUCHI_LOWER_MASS_LIMIT) >= 0) && (utils::Compare(m_Star2->MZAMS(), SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) {
+        double coreMass1 = m_Star1->MainSequenceCoreMass();
+        double coreMass2 = m_Star2->MainSequenceCoreMass();
+
+        double coreHeliumMass1 = m_Star1->HeliumAbundanceCore() * coreMass1;
+        double coreHeliumMass2 = m_Star2->HeliumAbundanceCore() * coreMass2;
+
+        finalHydrogenMass = finalMass * initialHydrogenFraction - coreHeliumMass1 - coreHeliumMass2;
+    }
+
     m_Star1->UpdateAfterMerger(finalMass, finalHydrogenMass);
 
     m_Star2->SwitchTo(STELLAR_TYPE::MASSLESS_REMNANT);
@@ -2078,7 +2088,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
             }
             else {                                                                                                              // have required mass loss
                 massDiffDonor = -massDiffDonor;                                                                                 // set mass difference
-                m_Donor->UpdateMinimumCoreMass();                                                                               // reset the minimum core mass following case A
+                m_Donor->UpdateTotalMassLossRate(massDiffDonor / (p_Dt * MYR_TO_YEAR));                                         // update mass loss rate for MS donor
+                m_Donor->UpdateMainSequenceCoreMass(p_Dt, massDiffDonor / (p_Dt * MYR_TO_YEAR));                                // update core mass for MS donor
             }
         }
 
@@ -2089,6 +2100,9 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
 
             double omegaDonor_pre_MT = m_Donor->Omega();                                                                        // used if full donor envelope is removed
 
+            m_Accretor->UpdateTotalMassLossRate(massGainAccretor / (p_Dt * MYR_TO_YEAR));                                       // update mass gain rate for MS accretor
+            m_Accretor->UpdateMainSequenceCoreMass(p_Dt, massGainAccretor / (p_Dt * MYR_TO_YEAR));                              // update core mass for MS accretor
+
             m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor);                                                 // set new mass of donor
             m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor);                                           // set new mass of accretor
 

src/BaseStar.cpp

@@ -129,12 +129,13 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
     m_Dt                                       = DEFAULT_INITIAL_DOUBLE_VALUE;
     m_Tau                                      = DEFAULT_INITIAL_DOUBLE_VALUE;
     m_Age                                      = 0.0;                                               // ensure age = 0.0 at construction (rather than default initial value)
+    m_MainSequenceCoreMass                     = DEFAULT_INITIAL_DOUBLE_VALUE;
     m_Mass                                     = m_MZAMS;
     m_Mass0                                    = m_MZAMS;
-    m_MinimumCoreMass                          = 0.0;
     m_Luminosity                               = m_LZAMS;
     m_Radius                                   = m_RZAMS;
     m_Temperature                              = m_TZAMS;
+    m_TotalMassLossRate                        = DEFAULT_INITIAL_DOUBLE_VALUE;
 	m_ComponentVelocity						   = Vector3d();
 
     m_OmegaCHE                                 = CalculateOmegaCHE(m_MZAMS, m_Metallicity);
@@ -2709,7 +2710,9 @@ double BaseStar::CalculateMassLossRate() {
 
         mDot = mDot * OPTIONS->OverallWindMassLossMultiplier();                                                     // apply overall wind mass loss multiplier
     }
-
+
+    UpdateTotalMassLossRate(-mDot);                                                                                 // update total mass loss rate
+
     return mDot;
 }
 
@@ -4571,12 +4574,15 @@ STELLAR_TYPE BaseStar::EvolveOnPhase(const double p_DeltaTime) {
     STELLAR_TYPE stellarType = m_StellarType;
 
     if (ShouldEvolveOnPhase()) {                                                    // evolve timestep on phase
+
+        UpdateMainSequenceCoreMass(p_DeltaTime, -m_Mdot);                           // update core mass, relevant for MS stars
+
         m_Tau        = CalculateTauOnPhase();
 
         m_COCoreMass = CalculateCOCoreMassOnPhase();
         m_CoreMass   = CalculateCoreMassOnPhase();
         m_HeCoreMass = CalculateHeCoreMassOnPhase();
-        
+
         m_Luminosity = CalculateLuminosityOnPhase();
 
         // Calculate abundances
@@ -4627,7 +4633,7 @@ STELLAR_TYPE BaseStar::ResolveEndOfPhase(const bool p_ResolveEnvelopeLoss) {
         if (p_ResolveEnvelopeLoss) stellarType = ResolveEnvelopeLoss();         // if required, resolve envelope loss if it occurs
 
         if (stellarType == m_StellarType) {                                     // staying on phase?
-
+            
             m_Tau         = CalculateTauAtPhaseEnd();
 
             m_COCoreMass  = CalculateCOCoreMassAtPhaseEnd();

src/BaseStar.h

@@ -137,13 +137,13 @@ class BaseStar {
             double              LogMetallicitySigma() const                                     { return m_Log10Metallicity; }                  // sigma in Hurley+ 2000
             double              LogMetallicityXi() const                                        { return m_Log10Metallicity - LOG10_ZSOL; }     // xi in Hurley+ 2000
             double              Luminosity() const                                              { return m_Luminosity; }
+            double              MainSequenceCoreMass() const                                    { return m_MainSequenceCoreMass; }
             double              Mass() const                                                    { return m_Mass; }
             double              Mass0() const                                                   { return m_Mass0; }
             double              MassPrev() const                                                { return m_MassPrev; }
             ST_VECTOR           MassTransferDonorHistory() const                                { return m_MassTransferDonorHistory; }
             double              Mdot() const                                                    { return m_Mdot; }
             double              Metallicity() const                                             { return m_Metallicity; }
-            double              MinimumCoreMass() const                                         { return m_MinimumCoreMass; }
             double              MZAMS() const                                                   { return m_MZAMS; }
             double              Omega() const                                                   { return m_AngularMomentum / CalculateMomentOfInertiaAU(); }
             double              OmegaCHE() const                                                { return m_OmegaCHE; }
@@ -185,6 +185,7 @@ class BaseStar {
             double              Temperature() const                                             { return m_Temperature; }
             double              Time() const                                                    { return m_Time; }
             double              Timescale(TIMESCALE p_Timescale) const                          { return m_Timescales[static_cast<int>(p_Timescale)]; }
+            double              TotalMassLossRate() const                                       { return m_TotalMassLossRate; }
             double              TZAMS() const                                                   { return m_TZAMS; }
     virtual ACCRETION_REGIME    WhiteDwarfAccretionRegime() const                               { return ACCRETION_REGIME::ZERO; }
             double              XExponent() const                                               { return m_XExponent; }
@@ -362,8 +363,9 @@ class BaseStar {
                                                        const double p_MassGainPerTimeStep,
                                                        const double p_Epsilon) { }                                                                                                  // Default is NO-OP
 
-    virtual void            UpdateMinimumCoreMass() { }                                                                                                                             // Only set minimal core mass following Main Sequence mass transfer to MS age fraction of TAMS core mass; default is NO-OP
-
+    virtual void            UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { }                                                                     // Set core mass for Main Sequence stars; default is NO-OP
+
+    virtual void            UpdateTotalMassLossRate(const double p_MassLossRate)                                { m_TotalMassLossRate = p_MassLossRate; }                           // m_TotalMassLossRate = -m_Mdot in SSE, during a mass transfer episode m_TotalMassLossRate = m_MassLossRateInRLOF
 
     // printing functions
     bool PrintDetailedOutput(const int p_Id, const SSE_DETAILED_RECORD_TYPE p_RecordType) const { 
@@ -434,9 +436,9 @@ class BaseStar {
     double                  m_HydrogenAbundanceSurface;                 // Hydrogen abundance at the surface
     bool                    m_LBVphaseFlag;                             // Flag to know if the star satisfied the conditions, at any point in its evolution, to be considered a Luminous Blue Variable (LBV)
     double                  m_Luminosity;                               // Current luminosity (Lsol)
+    double                  m_MainSequenceCoreMass;                     // Core mass of main sequence stars (Msol)
     double                  m_Mass;                                     // Current mass (Msol)
     double                  m_Mass0;                                    // Current effective initial mass (Msol)
-    double                  m_MinimumCoreMass;                          // Minimum core mass at end of main sequence (MS stars have no core in the Hurley prescription)
     double                  m_MinimumLuminosityOnPhase;                 // Only required for CHeB stars, but only needs to be calculated once per star
     double                  m_Mdot;                                     // Current mass loss rate in winds (Msol per yr)
     MASS_LOSS_TYPE          m_DominantMassLossRate;                     // Current dominant type of wind mass loss
@@ -446,6 +448,7 @@ class BaseStar {
     double                  m_Tau;                                      // Relative time
     double                  m_Temperature;                              // Current temperature (Tsol)
     double                  m_Time;                                     // Current physical time the star has been evolved (Myr)
+    double                  m_TotalMassLossRate;                        // Current mass loss/gain rate from mass transfer or winds (Msol per yr)
 
     // Previous timestep variables
     double                  m_DtPrev;                                   // Previous timestep

src/CH.h

@@ -88,7 +88,7 @@ class CH: virtual public BaseStar, public MS_gt_07 {
 
     bool            ShouldEvolveOnPhase() const                         { return m_Age < m_Timescales[static_cast<int>(TIMESCALE::tMS)] && (OPTIONS->OptimisticCHE() || Omega() >= m_OmegaCHE); } // Evolve on CHE phase if age in MS timescale and spinning at least as fast as CHE threshold
 
-    void            UpdateAgeAfterMassLoss();  
+    void            UpdateAgeAfterMassLoss();
 
 };
 

src/GiantBranch.h

@@ -135,7 +135,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
 
             void            UpdateInitialMass() { }                                                                                                                             // NO-OP for most stellar types
 
-            void            UpdateMinimumCoreMass() { }                                                                                                                         // NO-OP for most stellar types
+            void            UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { }                                                                 // NO-OP for most stellar types
 
 };
 

src/HG.cpp

@@ -818,6 +818,11 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const
 #define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)]      // for convenience and readability - undefined at end of function
 
     double RTMS = MainSequence::CalculateRadiusAtPhaseEnd(p_Mass, p_RZAMS);
+
+    // This ensures continuity of stellar tracks when SHIKAUCHI core mass prescription is used
+    if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0))
+        RTMS = MainSequence::CalculateRadiusAtPhaseEnd(m_Mass, p_RZAMS);
+
     double RGB  = GiantBranch::CalculateRadiusOnPhase_Static(p_Mass, m_Luminosity, b);
 
     double rx   = RGB;                                                                                              // Hurley sse terminlogy (rx)

src/HG.h

@@ -53,9 +53,9 @@ class HG: virtual public BaseStar, public GiantBranch {
         m_Age = m_Timescales[static_cast<int>(TIMESCALE::tMS)];                                                                                                                 // Set age appropriately
 
         // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass
-        // (only relevant if RetainCoreMassDuringCaseAMassTransfer())
-        if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) {
-            double desiredCoreMass = std::min(m_Mass, MinimumCoreMass());                                                                                                       // desired core mass
+        // (only relevant if MANDEL or SHIKAUCHI main sequence core mass prescription is used)
+        if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MainSequenceCoreMass())) < 0) {
+            double desiredCoreMass = std::min(m_Mass, MainSequenceCoreMass());                                                                                                       // desired core mass
             m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass);                                                                                                       // use root finder to find new core mass estimate
             if (m_Mass0 <= 0.0) {                                                                                                                                               // no root found - no solution for estimated core mass
                 m_Mass0 = m_Mass;                                                                                                                                               // if no root found we keep m_Mass0 equal to the total mass