First implementation of HAMSTARS mass transfer efficiency

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Tue Aug 19 15:39:23 2025 by COMPAS v03.25.00
+##~!!~## File Created Tue Sep  2 22:42:52 2025 by COMPAS v03.26.00
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -264,7 +264,7 @@ stringChoices:
 #    --critical-mass-ratio-prescription: 'NONE'                            # Default: 'NONE'                        # Options: ['HURLEY_HJELLMING_WEBBINK','GE','GE_IC','CLAEYS','NONE']
 #    --stellar-zeta-prescription: 'SOBERMAN'                               # Default: 'SOBERMAN'                    # Options: ['ARBITRARY','HURLEY','SOBERMAN']
 #    --mass-transfer-angular-momentum-loss-prescription: 'ISOTROPIC'       # Default: 'ISOTROPIC'                   # Options: ['ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS']
-#    --mass-transfer-accretion-efficiency-prescription: 'THERMAL'          # Default: 'THERMAL'                     # Options: ['FIXED','THERMAL']
+#    --mass-transfer-accretion-efficiency-prescription: 'THERMAL'          # Default: 'THERMAL'                     # Options: ['HAMSTARS','FIXED','THERMAL']
 #    --mass-transfer-rejuvenation-prescription: 'STARTRACK'                # Default: 'STARTRACK'                   # Options: ['STARTRACK','HURLEY']
 #    --mass-transfer-thermal-limit-accretor-multiplier: 'CFACTOR'          # Default: 'CFACTOR'                     # Options: ['ROCHELOBE','CFACTOR']
 #    --response-to-spin-up: 'TRANSFER_TO_ORBIT'                            # Default: 'TRANSFER_TO_ORBIT'           # Options: ['NO_LIMIT','KEPLERIAN_LIMIT','TRANSFER_TO_ORBIT']

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

@@ -887,7 +887,7 @@ DEPRECATION NOTICE: this option has been deprecated and will soon be removed. Pl
 
 **--mass-transfer-accretion-efficiency-prescription** |br|
 Mass transfer accretion efficiency prescription. |br|
-Options: { THERMAL, FIXED } |br|
+Options: { THERMAL, FIXED, HAMSTARS } |br|
 Default = THERMAL
 
 **--mass-transfer-angular-momentum-loss-prescription** |br|

online-docs/pages/whats-new.rst

@@ -3,6 +3,10 @@ What's new
 
 Following is a brief list of important updates to the COMPAS code.  A complete record of changes can be found in the file ``changelog.h``.
 
+**03.26.00 September 2, 2025**
+
+* Added HAMSTARS mass transfer efficiency prescription
+
 **03.25.00 August 19, 2025**
 
 * Added KLENCKI_LINEAR AM loss, which is linear in the specific AM gamma instead of the orbital separation (as in MACLEOD_LINEAR).

src/BaseBinaryStar.cpp

@@ -2092,7 +2092,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
     // can the mass transfer happen on a nuclear timescale?
     if (m_Donor->IsOneOf(NON_COMPACT_OBJECTS)) {
         // if MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION, then CalculateMassAcceptanceRate() already computed the correct m_FractionAccreted and massDiffDonor (no difference between nuclear and thermal timescale MT)
-        if (OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED) {
+        if (OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED || OPTIONS->MassTransferAccretionEfficiencyPrescription() == MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS) {
             // technically, we do not know how much mass the accretor should gain until we do the calculation,
             // which impacts the RL size, so we will check whether a nuclear timescale MT was feasible later
             massDiffDonor = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, -1.0, m_Dt);            // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, estimating accretion efficiency based on a mass donation rate of massDiffDonor/m_Dt for self-consistency

src/BaseStar.cpp

@@ -2713,12 +2713,17 @@ DBL_DBL BaseStar::CalculateMassAcceptanceRate(const double p_DonorMassRate, cons
 
     switch (OPTIONS->MassTransferAccretionEfficiencyPrescription()) {
 
-        case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED:                                   // thermally limited mass transfer:
-
+        case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED:                                   // thermally limited mass transfer
             acceptanceRate   = min(OPTIONS->MassTransferCParameter() * p_AccretorMassRate, p_DonorMassRate);
             fractionAccreted = acceptanceRate / p_DonorMassRate;
             break;
 
+        case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS:                                            // thermally limited mass transfer, following Lau+, 2024
+            // the mass transfer C parameter is fit using the data from Figure (1) of Lau+, 2024
+            acceptanceRate   = min(PPOW(10.0, (4.0 / PPOW((Mass() + 0.2), 0.3) - 0.6)) * p_AccretorMassRate, p_DonorMassRate);
+            fractionAccreted = acceptanceRate / p_DonorMassRate;
+            break;
+
         case MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION:                                      // fixed fraction of mass accreted, as in StarTrack
             fractionAccreted = OPTIONS->MassTransferFractionAccreted();
             acceptanceRate   = min(p_DonorMassRate, fractionAccreted * p_DonorMassRate);

src/CH.cpp

@@ -286,36 +286,33 @@ void CH::UpdateAgeAfterMassLoss() {
 
 
 /*
- * Calculate the weight for the OB and WR mass loss prescriptions
+ * CalculateMassLossFractionOB
  *
- * According to the prescription described in Yoon et al. 2006 (and Szecsi et al. 2015)
- * Use OB mass loss rates until Y1 = 0.55, WR mass loss rates above Y2 = 0.7, and 
- * linearly interpolate for Y1 < Ys < Y2 where Ys is the surface helium fraction.
- * 
- * Since we don't have Ys by default in our models, and detailed models show Ys
- * rises ~linearly from 0.2 to 1.0 across the main-sequence, here we simply 
- * use tau = t / tMS as a proxy for Ys.
+ * @brief
+ * Calculate the fraction of mass loss attributable to OB mass loss, per Yoon et al. 2006
+ *
+ * The model described in Yoon et al. 2006 (also Szecsi et al. 2015) uses OB mass loss while the
+ * He surface abundance is below 0.55, WR mass loss when the surface He abundance is above 0.7,
+ * and linearly interpolate when the He surface abundance is between those limits.
+ *
+ * This function calculates the fraction of mass loss attributable to OB mass loss, based on
+ * the He surface abundance and the abundance limits described in Yoon et al. 2006.  The value
+ * returned will be 1.0 if 100% of the mass loss is attributable to OB mass lass, 0.0 if 100% of
+ * the mass loss is attributable to WR mass loss, and in the range (0.0, 1.0) if the mass loss is
+ * a mix of OB and WR.
  *
- * CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface)
  *
- * @param   [IN]        p_HeliumAbundanceSurface           Surface helium abundance Ys
+ * double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const
  *
- * @return                                   weight for OB mass loss rate
+ * @param       p_HeAbundanceSurface            Helium abundance at the surface of the star
+ * @return                                      Fraction of mass loss attributable to OB mass loss
  */
-double CH::CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface) {
-
-    const double y1 = 0.55;
-    const double y2 = 0.70;
-    double       Ys = p_HeliumAbundanceSurface; // Get surface helium abundance
-    double weight   = 0.0;
-
+double CH::CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const {
 
-    // Calculate the weight as a function of Ys according to the prescription from Yoon et al. 2006
-         if (Ys < y1) weight = 1.0;
-    else if (Ys > y2) weight = 0.0;
-    else              weight = (y2 - Ys) / y2 - y1;
+    constexpr double limOB = 0.55;                                          // per Yoon et al. 2006
+    constexpr double limWR = 0.70;                                          // per Yoon et al. 2006
 
-    return weight;
+    return std::min(1.0, std::max (0.0, (limWR - p_HeAbundanceSurface) / (limWR - limOB)));
 }
 
 
@@ -337,7 +334,7 @@ double CH::CalculateMassLossRateBelczynski2010() {
     double Mdot   = 0.0;
     double MdotOB = 0.0;
     double MdotWR = 0.0;
-    double weight = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default 
+    double fractionOB = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default
 
     // Calculate OB mass loss rate according to the Vink et al. formalism
     MdotOB = BaseStar::CalculateMassLossRateBelczynski2010();
@@ -348,16 +345,16 @@ double CH::CalculateMassLossRateBelczynski2010() {
         // Calculate WR mass loss rate
         MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0);
 
-        // Calculate weight for combining these into total mass-loss rate
-        weight = CalculateMassLossRateWeightOB(m_HeliumAbundanceSurface);
+        // Calculate fraction for combining these into total mass-loss rate
+        fractionOB = CalculateMassLossFractionOB(m_HeliumAbundanceSurface);
 
     }
 
+    // Finally, combine each of these prescriptions according to the OB wind fraction
+    Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR);
+
     // Set dominant mass loss rate
-    m_DominantMassLossRate = weight == 0 ? MASS_LOSS_TYPE::WR : MASS_LOSS_TYPE::OB;   
-
-    // Finally, combine each of these prescriptions according to the weight
-    Mdot = (weight * MdotOB) + ((1.0 - weight) * MdotWR);
+    m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
 
     // Enhance mass loss rate due to rotation
     Mdot *= CalculateMassLossRateEnhancementRotation();
@@ -384,7 +381,7 @@ double CH::CalculateMassLossRateMerritt2025() {
     double Mdot   = 0.0;
     double MdotOB = 0.0;
     double MdotWR = 0.0;
-    double weight = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default
+    double fractionOB = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default
 
     // Calculate OB mass loss rate according to the chosen prescription
     MdotOB = BaseStar::CalculateMassLossRateOB(OPTIONS->OBMassLossPrescription());  
@@ -400,14 +397,14 @@ double CH::CalculateMassLossRateMerritt2025() {
         delete clone; clone = nullptr;                                                                          // return the memory allocated for the clone  
 
         // Calculate weight for combining these into total mass-loss rate
-        weight = CalculateMassLossRateWeightOB(m_HeliumAbundanceSurface);
+        fractionOB = CalculateMassLossFractionOB(m_HeliumAbundanceSurface);
     }
 
+    // Finally, combine each of these prescriptions according to the OB wind fraction
+    Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR);
+
     // Set dominant mass loss rate
-    m_DominantMassLossRate = weight == 0 ? MASS_LOSS_TYPE::WR : MASS_LOSS_TYPE::OB;   
-
-    // Finally, combine each of these prescriptions according to the weight
-    Mdot = (weight * MdotOB) + ((1.0 - weight) * MdotWR);
+    m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
 
     // Enhance mass loss rate due to rotation
     Mdot *= CalculateMassLossRateEnhancementRotation();

src/CH.h

@@ -74,7 +74,7 @@ class CH: virtual public BaseStar, public MS_gt_07 {
     // Mass loss rate
     double          CalculateMassLossRateBelczynski2010();
     double          CalculateMassLossRateMerritt2025();
-    double          CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface);
+    double          CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const;
 
     // Radius
     double          CalculateRadiusOnPhase() const                      { return m_RZAMS; }                                                                                         // Constant from birth

src/Options.cpp

@@ -2442,7 +2442,7 @@ std::string Options::OptionValues::CheckAndSetOptions() {
 
         if (m_UseMassTransfer && !DEFAULTED("mass-transfer-accretion-efficiency-prescription")) {                                   // mass transfer accretion efficiency prescription
             std::tie(found, m_MassTransferAccretionEfficiencyPrescription.type) = utils::GetMapKey(m_MassTransferAccretionEfficiencyPrescription.typeString, MT_ACCRETION_EFFICIENCY_PRESCRIPTION_LABEL, m_MassTransferAccretionEfficiencyPrescription.type);
-            COMPLAIN_IF(!found, "Unknown Mass Transfer Angular Momentum Loss Prescription");
+            COMPLAIN_IF(!found, "Unknown Mass Transfer Efficiency Prescription");
         }
 
         if (m_UseMassTransfer && !DEFAULTED("mass-transfer-angular-momentum-loss-prescription")) {                                  // mass transfer angular momentum loss prescription

src/changelog.h

@@ -1664,6 +1664,9 @@
 //                                          - Corrected calculation for Hurley Gamma constant C (C_GAMMA - see Hurley et al. 2000, just after eq 23, should use a(75) <= 1.0, not a(75) == 1.0 - confirmed in BSE Fortran source)
 //                                          - Added abs() to gamma calculation in Mainsequence.cpp::CalculateGamma() (per BSE Fortran source)
 //                                          - Clamped gamma to [0.0, gamma] in Mainsequence.cpp::CalculateGamma() (per discussion just after eq 23 - confirmed in BSE Fortran source)
+//  03.26.00  IM - September 2, 2025    - Enhancement, defect repairs:
+//                                          - First (simplified) implementation of the Lau+ (2024) Hamstars thermally limited accretion prescription
+//                                          - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010(), CalculateMassLossRateMerritt2025() and CH::CalculateMassLossFractionOB() [previously CalculateMassLossRateWeightOB()]
 //
 // Version string format is MM.mm.rr, where
 //
@@ -1674,7 +1677,7 @@
 // if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made
 // if mm is incremented, set rr to 00, even if defect repairs were also made
 
-const std::string VERSION_STRING = "03.25.02";
+const std::string VERSION_STRING = "03.26.00";
 
 
 # endif // __changelog_h__

src/typedefs.h

@@ -643,10 +643,11 @@ const COMPASUnorderedMap<METALLICITY_DISTRIBUTION, std::string> METALLICITY_DIST
 };
 
 // mass transfer accretion efficiency prescriptions
-enum class MT_ACCRETION_EFFICIENCY_PRESCRIPTION: int { THERMALLY_LIMITED, FIXED_FRACTION };
+enum class MT_ACCRETION_EFFICIENCY_PRESCRIPTION: int { THERMALLY_LIMITED, FIXED_FRACTION, HAMSTARS };
 const COMPASUnorderedMap<MT_ACCRETION_EFFICIENCY_PRESCRIPTION, std::string> MT_ACCRETION_EFFICIENCY_PRESCRIPTION_LABEL = {
     { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::THERMALLY_LIMITED, "THERMAL" },
-    { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION,    "FIXED" }
+    { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION,    "FIXED" },
+    { MT_ACCRETION_EFFICIENCY_PRESCRIPTION::HAMSTARS,          "HAMSTARS"}
 };
 
 // mass transfer angular momentum loss prescriptions