Updating mass transfer physics

Author: Ilya Mandel

online-docs/pages/whats-new.rst

@@ -4,8 +4,9 @@ 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.21.00 July 6, 2025**
-* - Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
+* Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
 * Added options to set beta and gamma prescription for second stage of 2-stage CE (``--common-envelope-second-stage-beta``, ``--common-envelope-second-stage-gamma-prescription``)
+* Fixed a bug in the calculation of zeta_equilibrium, which impacts when mass transfer is declared to proceed on a nuclear timescale (and hence how conservative it is)
 
 **03.20.06 June 25, 2025**
 

src/BaseBinaryStar.cpp

@@ -2058,27 +2058,22 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
     }
 
     m_MassTransferTimescale         = MT_TIMESCALE::NONE;                                                                       // initial reset
-    double betaThermal              = 0.0;                                                                                      // fraction of mass accreted if accretion proceeds on thermal timescale
-    double maximumAccretionRate     = 0.0;                                                                                      // accretion rate if accretion proceeds on thermal timescale
+
     double donorMassLossRateThermal = m_Donor->CalculateThermalMassLossRate();
 
-    std::tie(maximumAccretionRate, betaThermal) = m_Accretor->CalculateMassAcceptanceRate(donorMassLossRateThermal,
-                                                                                          m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius),
-                                                                                          donorIsHeRich);
-    double massDiffDonor = 0.0;
-        
+    std::tie(std::ignore, m_FractionAccreted) = m_Accretor->CalculateMassAcceptanceRate(donorMassLossRateThermal,
+                                                                                        m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), donorIsHeRich);
+    double massDiffDonor            = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted, 0.0);         // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, fixed beta
+            
     // can the mass transfer happen on a nuclear timescale?
     if (m_Donor->IsOneOf(NON_COMPACT_OBJECTS)) {
-        // 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
-        double maximumAccretedMass = maximumAccretionRate * m_Dt;
+        // 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) {
-            massDiffDonor      = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, -1.0, maximumAccretedMass);            // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, fixed accretion amount
-            m_FractionAccreted = std::min(maximumAccretedMass, massDiffDonor) / massDiffDonor;
-        }
-        else {
-            m_FractionAccreted = maximumAccretionRate / donorMassLossRateThermal;   // relevant for MT_ACCRETION_EFFICIENCY_PRESCRIPTION::FIXED_FRACTION
-            massDiffDonor      = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, m_FractionAccreted, 0.0);              // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe, fixed beta
+            // 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
+            std::tie(std::ignore, m_FractionAccreted) = m_Accretor->CalculateMassAcceptanceRate(massDiffDonor/m_Dt,
+                                                                                                m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), donorIsHeRich);
         }
 
         // check that the star really would have consistently fit into the Roche lobe
@@ -2091,13 +2086,11 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
             m_ZetaLobe              = zetaLobe;
         }
     }
-    if (m_MassTransferTimescale != MT_TIMESCALE::NUCLEAR) {                                                                     // thermal timescale mass transfer (we will check for dynamically unstable / CE mass transfer later)
-        m_ZetaLobe              = CalculateZetaRocheLobe(jLoss, betaThermal);
+    if (m_MassTransferTimescale != MT_TIMESCALE::NUCLEAR) {                                                                     // thermal timescale mass transfer (we will check for dynamically unstable / CE mass transfer later); m_FractionAccreted and massDiffDonor already computed for this case
+        m_ZetaLobe              = CalculateZetaRocheLobe(jLoss, m_FractionAccreted);
         m_ZetaStar              = m_Donor->CalculateZetaAdiabatic();
         m_MassLossRateInRLOF    = donorMassLossRateThermal;
-        m_FractionAccreted      = betaThermal;
         m_MassTransferTimescale = MT_TIMESCALE::THERMAL;
-        massDiffDonor           = MassLossToFitInsideRocheLobe(this, m_Donor, m_Accretor, betaThermal, 0.0);                    // use root solver to determine how much mass should be lost from the donor to allow it to fit within the Roche lobe
     }
 
     double aInitial = m_SemiMajorAxis;                                                                                          // semi-major axis in default units, AU, current timestep

src/BaseBinaryStar.h

@@ -577,13 +577,13 @@ class BaseBinaryStar {
      *
      *
      * Constructor: initialise the class
-     * template <class T> RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass, ERROR *p_Error)
+     * template <class T> RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_Dt, ERROR *p_Error)
      *
      * @param   [IN]    p_Binary                    (Pointer to) The binary star under examination
      * @param   [IN]    p_Donor                     (Pointer to) The star donating mass
      * @param   [IN]    p_Accretor                  (Pointer to) The star accreting mass
-     * @param   [IN]    p_FractionAccreted          The fraction of the donated mass accreted by the accretor (for thermal timescale accretion)
-     * @param   [IN]    p_MaximumAccretedMass       The total amount of mass that can be accreted (for nuclear timescale accretion, p_FractionAccreted should be negative for this to be used)
+     * @param   [IN]    p_FractionAccreted          The fraction of the donated mass accreted by the accretor (if known in advance, otherwise zero)
+     * @param   [IN]    p_Dt                        Time step duration (relevant for nuclear timescale mass transfer)
      * @param   [IN]    p_Error                     (Address of variable to record) Error encountered in functor
      * 
      * Function: calculate radius difference after mass loss
@@ -594,13 +594,13 @@ class BaseBinaryStar {
      */    
     template <class T>
     struct RadiusEqualsRocheLobeFunctor {
-        RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_MaximumAccretedMass, ERROR *p_Error) {
+        RadiusEqualsRocheLobeFunctor(BaseBinaryStar *p_Binary, BinaryConstituentStar *p_Donor, BinaryConstituentStar *p_Accretor, double p_FractionAccreted, double p_Dt, ERROR *p_Error) {
             m_Binary           = p_Binary;
             m_Donor            = p_Donor;
             m_Accretor         = p_Accretor;
             m_Error            = p_Error;
             m_FractionAccreted = p_FractionAccreted;
-            m_MaximumAccretedMass = p_MaximumAccretedMass;
+            m_Dt               = p_Dt;
         }
         T operator()(double const& p_dM) {
 
@@ -611,12 +611,16 @@ class BaseBinaryStar {
 
             double donorMass     = m_Donor->Mass();
             double accretorMass  = m_Accretor->Mass();
+            // use stale value of accretor RL radius -- this is only relevant for nuclear timescale MT, when the change in accretor RL radius should be small
+            double accretorRLradius = CalculateRocheLobeRadius_Static(accretorMass, donorMass) * AU_TO_RSOL * m_Binary->SemiMajorAxis() * (1.0 - m_Binary->Eccentricity());
 
             // beta is the actual accretion efficiency; if p_FractionAccreted is negative (placeholder
             // for nuclear timescale accretion efficiency, for which the total accretion mass over the
-            // duration of the timestep is known), then the ratio of the maximum allowed accreted
-            // mass / donated mass is used
-            double beta = (utils::Compare(m_FractionAccreted, 0.0) >=0 ) ? m_FractionAccreted : std::min(m_MaximumAccretedMass/p_dM, 1.0);
+            // duration of the timestep is known), then must estimate it on the fly for consistency
+            double beta = m_FractionAccreted;
+            if(utils::Compare(beta, 0.0) < 0)
+                std::tie(std::ignore, beta) = m_Accretor->CalculateMassAcceptanceRate(p_dM / m_Dt,
+                                                    m_Accretor->CalculateThermalMassAcceptanceRate(accretorRLradius), m_Donor->IsOneOf(He_RICH_TYPES));
 
             double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorMass, -p_dM , *m_Accretor, beta, false);
             double RLRadius      = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - p_dM, accretorMass + (beta * p_dM)) * AU_TO_RSOL;
@@ -631,7 +635,7 @@ class BaseBinaryStar {
         BinaryConstituentStar *m_Accretor;
         ERROR                 *m_Error;
         double                 m_FractionAccreted;
-        double                 m_MaximumAccretedMass;
+        double                 m_Dt;
     };
 
 

src/BaseStar.h

@@ -302,7 +302,7 @@ class BaseStar {
             double          CalculateZetaAdiabatic();
     virtual double          CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription)              { return 0.0; }                                                     // Use inheritance hierarchy
 
-            double          CalculateZetaEquilibrium()                                                          { return 0.0; }
+    virtual double          CalculateZetaEquilibrium()                                                          { return 0.0; }
 
             void            ClearCurrentSNEvent()                                                               { m_SupernovaDetails.events.current = SN_EVENT::NONE; }             // Clear supernova event/state for current timestep
 

src/changelog.h

@@ -1599,10 +1599,11 @@
 //                                          - Fixed error in MainSequence::CalculateInitialMainSequenceCoreMass()
 //  03.20.09  RTW - Jun 30, 2025        - Enhancement:
 //                                          - Added individual velocity components for stars to the LogTypedefs file so they can be included in the output (as ANY_STAR_PROPERTY::VELOCITY_X, or Y, Z)
-//  03.21.00   IM - July 6, 2025        - Enhancements:
+//  03.21.00   IM - July 6, 2025        - Enhancements, defect repair:
 //                                          - Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
 //                                          - Added options to set beta and gamma prescription for second stage of 2-stage CE (--common-envelope-second-stage-beta, --common-envelope-second-stage-gamma-prescription)
-
+//                                          - Fixed a bug in CalculateZetaEquilibrium(), which impacted when mass transfer is declared nuclear (and how conservative it is)
+//                                          - Now calculate mass accretion rate for nuclear timescale mass transfer on the fly to match with donor mass loss rate set by donor mass loss (required to fit into Roche lobe) divided by time step 
 
 
 const std::string VERSION_STRING = "03.21.00";