Merge pull request #1401 from brcekadam/CHE_and_merger_fixes

brcekadam · web-flow · commit 3ad655bf9c00 · 2025-06-23T16:40:46.000+10:00
Fixes to CHE and MS mergers with BRCEK core mass prescription
diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp
@@ -1778,15 +1778,14 @@ void BaseBinaryStar::ResolveMainSequenceMerger() {
         (utils::Compare(m_Star1->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)                &&
         (utils::Compare(m_Star2->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)) {
         
-        double coreMass1       = m_Star1->MainSequenceCoreMass();
-        double coreMass2       = m_Star2->MainSequenceCoreMass();
+        // total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope)
+        double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + m_Star1->HydrogenAbundanceSurface() * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0);
+        double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + m_Star2->HydrogenAbundanceSurface() * std::max(m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass(), 0.0);
         
-        double coreHeliumMass1 = m_Star1->HeliumAbundanceCore() * coreMass1;
-        double coreHeliumMass2 = m_Star2->HeliumAbundanceCore() * coreMass2;
-        
-        finalHydrogenMass      = finalMass * initialHydrogenFraction - coreHeliumMass1 - coreHeliumMass2;
+        // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope of the star with higher surface hydrogen abundance
+        finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * std::max(m_Star1->HydrogenAbundanceSurface(), m_Star2->HydrogenAbundanceSurface());
     }
-    else finalHydrogenMass     = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction;
+    else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction;
        
     m_SemiMajorAxis = std::numeric_limits<float>::infinity();                                   // set separation to infinity to avoid subsequent fake interactions with a massless companion (RLOF, CE, etc.)
     
diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp
@@ -107,6 +107,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
     // initialise remaining member variables
 
     // Zero age main sequence parameters
+    m_InitialMainSequenceCoreMass              = DEFAULT_INITIAL_DOUBLE_VALUE;                      // initialised in MS_gt_07 class if BRCEK core mass prescription is used
     m_RZAMS                                    = CalculateRadiusAtZAMS(m_MZAMS);
     m_LZAMS                                    = CalculateLuminosityAtZAMS(m_MZAMS);
     m_TZAMS                                    = CalculateTemperatureOnPhase_Static(m_LZAMS, m_RZAMS);
diff --git a/src/BaseStar.h b/src/BaseStar.h
@@ -106,6 +106,7 @@ class BaseStar {
     double              HydrogenAbundanceSurface() const                                { return m_HydrogenAbundanceSurface; }
     double              InitialHeliumAbundance() const                                  { return m_InitialHeliumAbundance; }
     double              InitialHydrogenAbundance() const                                { return m_InitialHydrogenAbundance; }
+    double              InitialMainSequenceCoreMass() const                             { return m_InitialMainSequenceCoreMass; }
     bool                IsAIC() const                                                   { return (m_SupernovaDetails.events.current & SN_EVENT::AIC) == SN_EVENT::AIC; }
     bool                IsCCSN() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::CCSN) == SN_EVENT::CCSN; }
     bool                IsHeSD() const                                                  { return (m_SupernovaDetails.events.current & SN_EVENT::HeSD) == SN_EVENT::HeSD; }
@@ -404,6 +405,7 @@ class BaseStar {
     // Zero Age Main Sequence
     double                  m_InitialHeliumAbundance;                   // Initial helium abundance (Y)
     double                  m_InitialHydrogenAbundance;                 // Initial hydrogen abundance (X)
+    double                  m_InitialMainSequenceCoreMass;              // Initial main sequence core mass (used in BRCEK core mass prescription)
     double                  m_LZAMS;                                    // ZAMS Luminosity
     double                  m_MZAMS;                                    // ZAMS Mass
     double                  m_OmegaZAMS;                                // ZAMS Angular Frequency
diff --git a/src/CH.cpp b/src/CH.cpp
@@ -423,6 +423,10 @@ STELLAR_TYPE CH::EvolveToNextPhase() {
     if (m_Age < m_Timescales[static_cast<int>(TIMESCALE::tMS)]) {           // evolving off because of age?
         stellarType = STELLAR_TYPE::MS_GT_07;                               // no - must have spun down - evolve as MS star now
         m_CHE       = false;                                                // evolved CH->MS
+        
+        // if BRCEK core mass calculations enabled, initialise the core mass based on current mass and central helium fraction
+        if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0))
+            m_MainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_Mass, m_HeliumAbundanceCore);
     }
     else {                                                                  // yes
         stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_MS;                   // evolve as HeMS star now
diff --git a/src/CH.h b/src/CH.h
@@ -89,6 +89,8 @@ 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            UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate)    { };                                                                         // Do not use core mass calculations during CHE phase
 
 };
 
diff --git a/src/MS_gt_07.h b/src/MS_gt_07.h
@@ -48,7 +48,7 @@ class MS_gt_07: virtual public BaseStar, public MainSequence {
             utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0                         &&                                                             // ZAMS mass >= BRCEK_LOWER_MASS_LIMIT?
             m_Time <= 0.0) {                                                                                                                            // star not yet aged past creation?
                                                                                                                                                         // yes - initialise
-            m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS);
+            m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS, m_InitialHeliumAbundance);
             m_MainSequenceCoreMass        = m_InitialMainSequenceCoreMass;
             m_Luminosity                  = MainSequence::CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0);
             m_Radius                      = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Tau, m_RZAMS0);
diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp
@@ -289,8 +289,9 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl
 #define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)]  // for convenience and readability - undefined at end of function
     
     // If BRCEK core prescription is used, return luminosity from Shikauchi et al. (2024) during core hydrogen burning (valid for MZAMS >= 15 Msol) or
-    // luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT)
-    if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) {
+    // luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT); do not use Shikauchi luminosity
+    // prescription during CHE
+    if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0) && !m_CHE) {
         if (utils::Compare(p_Time, 0.99 * timescales(tMS)) > 0)                                                                 // star in MS hook?
             return CalculateLuminosityTransitionToHG(p_Mass, p_Time, p_LZAMS);
         else {
@@ -826,18 +827,30 @@ DBL_DBL MainSequence::CalculateMainSequenceCoreMassBrcek(const double p_Dt, cons
 
 
 /*
- * Calculate the initial convective core mass of a main sequence star using Equation (A3) from Shikauchi et al. (2024),
- * also used for calculating core mass after MS merger
+ * Calculate the initial convective core mass of a main sequence star at ZAMS using Equation (A3) from Shikauchi+ (2024),
+ * or after full mixing (due to merger or CHE) for an arbitrary central helium fraction using the approach
+ * described in Brcek+ (2025)
  *
- * double CalculateInitialMainSequenceCoreMass(const double p_MZAMS)
+ * double CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore)
  *
- * @param   [IN]    p_MZAMS                     Mass at ZAMS or after merger in Msol
+ * @param   [IN]    p_Mass                      Mass at ZAMS, after merger or after spin down of CH star in Msol
+ * @param   [IN]    p_HeliumAbundanceCore       Central helium fraction
  * @return                                      Mass of the convective core at ZAMS or after merger in Msol
  */
-double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const {
-    DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS);
-    double fmix                 = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_MZAMS / fmixCoefficients[2]);
-    return fmix * p_MZAMS;
+double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) const {
+    
+    double fmix = 0.0;
+    // At ZAMS, use the approach from Shikauchi+ (2024)
+    if (utils::Compare(p_HeliumAbundanceCore, m_InitialHeliumAbundance) == 0) {
+        DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS);
+        fmix                        = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_Mass / fmixCoefficients[2]);
+    }
+    // After full mixing not at ZAMS, use the approach from Brcek+ (2025)
+    else {
+        double h = PPOW(10.0, p_HeliumAbundanceCore * (p_HeliumAbundanceCore + 2.0) / 4.0);
+        fmix     = BRCEK_FMIX_COEFFICIENTS[0] + BRCEK_FMIX_COEFFICIENTS[1] * std::exp(-p_Mass * h / BRCEK_FMIX_COEFFICIENTS[2]) * PPOW(1.0 - BRCEK_FMIX_COEFFICIENTS[4] / (p_Mass * h), BRCEK_FMIX_COEFFICIENTS[3]);
+    }
+    return fmix * p_Mass;
 }
 
 
@@ -1142,12 +1155,14 @@ void MainSequence::UpdateAfterMerger(double p_Mass, double p_HydrogenMass) {
     m_Age = m_Tau * timescales(tMS);
     
     m_HeliumAbundanceCore   = 1.0 - m_Metallicity - p_HydrogenMass / p_Mass;
-    
     m_HydrogenAbundanceCore = 1.0 - m_Metallicity - m_HeliumAbundanceCore;
     
+    m_HeliumAbundanceSurface   = m_HeliumAbundanceCore;                                         // abundances are the same throughout the star, assuming uniform mixing after merger
+    m_HydrogenAbundanceSurface = m_HydrogenAbundanceCore;
+    
     if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) {
-        m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass);           // update initial mixing core mass
-        m_MainSequenceCoreMass        = m_InitialMainSequenceCoreMass;                          // update core mass
+        m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass, m_HeliumAbundanceCore);           // update initial mixing core mass
+        m_MainSequenceCoreMass        = m_InitialMainSequenceCoreMass;                                                 // update core mass
     }
     
     UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true);
diff --git a/src/MainSequence.h b/src/MainSequence.h
@@ -25,9 +25,7 @@ class MainSequence: virtual public BaseStar {
 protected:
     
     // member variables
-    
     double          m_HeliumAbundanceCoreOut      = m_InitialHeliumAbundance;                                                                                       // Helium abundance just outside the core, used for rejuvenation calculations
-    double          m_InitialMainSequenceCoreMass = 0.0;                                                                                                            // Initial mass of the mixing core is initialised in MS_gt_07 class
 
     // member functions - alphabetically
     double          CalculateAlphaL(const double p_Mass) const;
@@ -78,7 +76,7 @@ class MainSequence: virtual public BaseStar {
     double          CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore) const;
     double          CalculateLuminosityTransitionToHG(const double p_Mass, const double p_Age, double const p_LZAMS) const;
     DBL_DBL         CalculateMainSequenceCoreMassBrcek(const double p_Dt, const double p_MassLossRate);
-    double          CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const;
+    double          CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) const;
     double          CalculateMomentOfInertia() const                                        { return (0.1 * (m_Mass) * m_Radius * m_Radius); }                      // k2 = 0.1 as defined in Hurley et al. 2000, after eq 109
 
     double          CalculatePerturbationMu() const                                         { return 5.0; }                                                         // mu(MS) = 5.0 (Hurley et al. 2000, eqs 97 & 98)
diff --git a/src/Star.h b/src/Star.h
@@ -97,6 +97,7 @@ class Star {
     double              HydrogenAbundanceSurface() const                                                            { return m_Star->HydrogenAbundanceSurface(); }
     double              InitialHeliumAbundance() const                                                              { return m_Star->InitialHeliumAbundance(); }
     double              InitialHydrogenAbundance() const                                                            { return m_Star->InitialHydrogenAbundance(); }
+    double              InitialMainSequenceCoreMass() const                                                         { return m_Star->InitialMainSequenceCoreMass(); }
     bool                IsAIC() const                                                                               { return m_Star->IsAIC(); }
     bool                IsCCSN() const                                                                              { return m_Star->IsCCSN(); }
     bool                IsDegenerate() const                                                                        { return m_Star->IsDegenerate(); }
diff --git a/src/changelog.h b/src/changelog.h
@@ -1583,8 +1583,11 @@
 //                                          - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02)
 //                                          - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static()
 //                                          - Removed McSN from GBParams, instead computed on the fly when needed
+//  03.20.04   AB - Jun 23, 2025        - Defect repair, enhancement:
+//                                          - Fixes to MS mergers and CHE when BRCEK core mass prescription is used -- MS core mass is now correctly initialised after full mixing in MS
+//                                            mergers and CH stars that spun down
 
-const std::string VERSION_STRING = "03.20.03";
+const std::string VERSION_STRING = "03.20.04";
 
 
 # endif // __changelog_h__
diff --git a/src/constants.h b/src/constants.h
@@ -3811,5 +3811,9 @@ const std::vector<DBL_VECTOR> SHIKAUCHI_L_COEFFICIENTS = {
     {3.45464814, 0.94880846, -1.11409154, -0.86672079, 2.38986855, 0.04448855, -2.74913945, 0.60905625, -0.27648361,  0.03514139, 1.37569819, -0.19184532, 0.12816567, -0.14392935, 1.76390159},   // 1/3*Z_Sun
     {3.80166901, 0.37407948, -1.29904749, -1.34541622, 3.70934166, 0.28320469, -3.92327169, 0.92444477, -0.40146717, -0.00821364, 1.80297947, -0.15776603, 0.09205681, -0.21913557, 1.78496679}    // Solar metallicity Z_Sun
 };
+// Coefficients used to determine the initial convective core mass of MS star after full mixing (due to a merger or CHE)
+// from Brcek et al. (2025)
+const DBL_VECTOR BRCEK_FMIX_COEFFICIENTS = {0.897487516925644, -0.570499075062176, 56.9817940103868, 0.54483109433479, 1.72842835324164};
+
 
 #endif // __constants_h__
