Updates to WD accretion to address issue #1384

online-docs/pages/whats-new.rst

@@ -3,13 +3,17 @@ 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.20.01 May 26, 2025**
+
+* Updates to mass accretion for massive ONe WD 
+* Changed white dwarf mass-radius relation to use expression from Eggleton 1986, suitable for extremely low-mass white dwarfs.
+
 **03.20.00 May 25, 2025**
 
 * Replaced the name of the ``KAPIL2024`` tides prescription with ``KAPIL2025``.
 * Updated the equilibrium and dynamical tides equations to match the paper.
 * New outputs for BSE_DETAILED_OUTPUT from tidal evolution, including ``CIRCULARIZATION_TIMESCALE``, ``SYNCHRONIZATION_TIMESCALE_1``, ``SYNCHRONIZATION_TIMESCALE_2``, ``TIDAL_POTENTIAL_LOVE_NUMBER_22_1``, ``TIDAL_POTENTIAL_LOVE_NUMBER_10_EQ_1``, and ``TIDAL_POTENTIAL_LOVE_NUMBER_32_DYN_2``
 
-
 **03.19.00 May 22, 2025**
 
 * Added functionality to create new System Snapshot logfile |br|

src/BaseBinaryStar.cpp

@@ -1673,8 +1673,8 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
             THROW_ERROR(ERROR::UNKNOWN_CE_FORMALISM);                                                                   // throw error
     }
 
-	double rRLdfin1     = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final);                // Roche-lobe radius in AU after CEE, seen by star1
-	double rRLdfin2     = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final);                // Roche-lobe radius in AU after CEE, seen by star2
+    double rRLdfin1     = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass1Final, m_Mass2Final);                // Roche-lobe radius in AU after CEE, seen by star1
+    double rRLdfin2     = m_SemiMajorAxis * CalculateRocheLobeRadius_Static(m_Mass2Final, m_Mass1Final);                // Roche-lobe radius in AU after CEE, seen by star2
     double rRLdfin1Rsol = rRLdfin1 * AU_TO_RSOL;                                                                        // Roche-lobe radius in Rsol after CEE, seen by star1
     double rRLdfin2Rsol = rRLdfin2 * AU_TO_RSOL;                                                                        // Roche-lobe radius in Rsol after CEE, seen by star2
     m_Eccentricity      = 0.0;                                                                                          // we assume that a common envelope event (CEE) circularises the binary
@@ -1684,13 +1684,13 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
 
     // update stellar type after losing its envelope. Star1, Star2 or both if double CEE.
 
-    if (isDonorMS || (!envelopeFlag1 && !envelopeFlag2)) {                                                              // stellar merger
+    if (isDonorMS || (!envelopeFlag1 && !envelopeFlag2) || (m_Star1->IsRLOF() && m_Star1->IsOneOf(WHITE_DWARFS)) || (m_Star2->IsRLOF() && m_Star2->IsOneOf(WHITE_DWARFS))) {                                                              // stellar merger
         m_MassTransferTrackerHistory = MT_TRACKING::MERGER; 
         m_Flags.stellarMerger        = true;
     }
     else if ((m_Star1->DetermineEnvelopeType()==ENVELOPE::RADIATIVE && !m_Star1->IsOneOf(ALL_MAIN_SEQUENCE)) ||
              (m_Star2->DetermineEnvelopeType()==ENVELOPE::RADIATIVE && !m_Star2->IsOneOf(ALL_MAIN_SEQUENCE)) ) {        // check if we have a non-MS radiative-envelope star
-        if (!OPTIONS->AllowRadiativeEnvelopeStarToSurviveCommonEnvelope() && OPTIONS->CommonEnvelopeFormalism()!=CE_FORMALISM::TWO_STAGE) { // stellar merger
+        if (!OPTIONS->AllowRadiativeEnvelopeStarToSurviveCommonEnvelope() && OPTIONS->CommonEnvelopeFormalism() != CE_FORMALISM::TWO_STAGE) { // stellar merger
             m_CEDetails.optimisticCE     = true;
             m_MassTransferTrackerHistory = MT_TRACKING::MERGER;
             m_Flags.stellarMerger        = true;
@@ -1702,7 +1702,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
         m_Flags.stellarMerger = true;
     }
 
-	if (!m_Flags.stellarMerger) {                                                                                       // stellar merger?
+    if (!m_Flags.stellarMerger) {                                                                                       // stellar merger?
                                                                                                                         // no - continue evolution
         STELLAR_TYPE stellarType1 = m_Star1->StellarType();                                                             // star 1 stellar type before resolving envelope loss
         STELLAR_TYPE stellarType2 = m_Star2->StellarType();                                                             // star 2 stellar type before resolving envelope loss
@@ -1731,7 +1731,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {
 
         m_Star1->SetPostCEEValues();                                                                                    // squirrel away post CEE stellar values for star 1
         m_Star2->SetPostCEEValues();                                                                                    // squirrel away post CEE stellar values for star 2
-        SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, semiMajorAxisAfterStage1 * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol);                     // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end)
+        SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, semiMajorAxisAfterStage1 * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end)
 
         if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) {// is there immediate post-CE RLOF which is not allowed?
             m_MassTransferTrackerHistory = MT_TRACKING::MERGER;

src/BaseStar.h

@@ -300,8 +300,8 @@ class BaseStar {
 
             void            ClearCurrentSNEvent()                                                               { m_SupernovaDetails.events.current = SN_EVENT::NONE; }             // Clear supernova event/state for current timestep
 
-    virtual ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich,
-                                                      const double p_DonorThermalMassLossRate)                  { return ACCRETION_REGIME::ZERO; }                                  // Placeholder, use inheritance for WDs
+    virtual ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate,
+                                                      const bool p_HeRich)                                      { return ACCRETION_REGIME::ZERO; }                                  // Placeholder, use inheritance for WDs
 
     virtual ENVELOPE        DetermineEnvelopeType() const                                                       { return ENVELOPE::REMNANT; }                                       // Default is REMNANT - but should never be called
     virtual MT_CASE         DetermineMassTransferTypeAsDonor() const                                            { return MT_CASE::OTHER; }                                          // Not A, B, C, or NONE

src/CHeB.h

@@ -76,7 +76,7 @@ class CHeB: virtual public BaseStar, public FGB {
     double          CalculateCoreMassOnPhase(const double p_Mass, const double p_Tau) const;
     double          CalculateCoreMassOnPhase() const                            { return CalculateCoreMassOnPhase(m_Mass0, m_Tau); }                            // Use class member variables
 
-    double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const    { return 0.33; }                                                                // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
+    double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const    { return HURLEY_HJELLMING_WEBBINK_QCRIT_MS_GT_07; }
 
     double          CalculateHeCoreMassAtPhaseEnd() const                       { return m_CoreMass; }
 

src/COWD.cpp

@@ -1,105 +1,5 @@
 #include "COWD.h"
 
-/* For COWDs, calculate:
- *
- *     (a) the maximum mass acceptance rate of this star, as the accretor, during mass transfer, and
- *     (b) the retention efficiency parameter
- *
- *
- * For a given mass transfer rate, this function computes the amount of mass a WD would retain after
- * flashes, as given by appendix B of Claeys+ 2014. 
- * https://ui.adsabs.harvard.edu/abs/2014A%26A...563A..83C/abstract 
- *
- *
- * DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich)
- *
- * @param   [IN]    p_DonorMassRate             Mass transfer rate from the donor
- * @param   [IN]    p_IsHeRich                  Material is He-rich or not
- * @return                                      Tuple containing the Maximum Mass Acceptance Rate (Msun/yr) and Retention Efficiency Parameter
- */
-DBL_DBL COWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) {
-
-    m_AccretionRegime = DetermineAccretionRegime(p_IsHeRich, p_DonorMassRate); 
-
-    double acceptanceRate   = 0.0;                                                       // acceptance mass rate - default = 0.0
-    double fractionAccreted = 0.0;                                                       // accretion fraction - default = 0.0
-
-    acceptanceRate = p_DonorMassRate * CalculateEtaHe(p_DonorMassRate);
-    if (!p_IsHeRich) acceptanceRate *= CalculateEtaH(p_DonorMassRate);
-
-    fractionAccreted = acceptanceRate / p_DonorMassRate;
-
-    return std::make_tuple(acceptanceRate, fractionAccreted);
-}
-
-
-/* 
- * Determine the WD accretion regime based on the MT rate and whether the donor is He rich. Also,
- * initialize He-Shell detonation or Off-center ignition when necessary, by changing the value
- * of m_HeShellDetonation or m_OffCenterIgnition (respectively).
- *
- * The accretion regime is one of the options listed in enum ACCRETION_REGIME (constants.h)
- *
- * Note that we have merged the different flashes regimes from Piersanti+ 2014 into a single regime.
- *
- * ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, const double p_DonorMassLossRate) 
- *
- * @param   [IN]    p_HeRich                 Whether the accreted material is helium-rich or not
- * @param   [IN]    p_DonorMassLossRate      Donor mass loss rate, in units of Msol / Myr
- * @return                                   Current WD accretion regime
- */
-ACCRETION_REGIME COWD::DetermineAccretionRegime(const bool p_HeRich, const double p_DonorMassLossRate) {
-
-    double logMdot          = log10(p_DonorMassLossRate / MYR_TO_YEAR);                                                     // logarithm of the accreted mass (M_sun/yr)
-    ACCRETION_REGIME regime = ACCRETION_REGIME::ZERO;
-
-    if (p_HeRich) {
-        // The following coefficients in logMassTransfer limits come from table A1 in Piersanti+ 2014.
-        double logMassTransferCrit       = WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_0 + WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_1 * m_Mass;
-        double logMassTransferStable     = WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_0 + WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_1 * m_Mass;  // Piersanti+2014 has several Flashes regimes. Here we group them into one.
-        double logMassTransferDetonation = WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_0 + WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_1 * m_Mass;  // critical value for double detonation regime in Piersanti+ 2014
-        if (utils::Compare(logMdot, logMassTransferStable) < 0) {
-            if (utils::Compare(logMdot, logMassTransferDetonation) > 0) {
-                regime = ACCRETION_REGIME::HELIUM_FLASHES;
-            } 
-            else {
-                regime = ACCRETION_REGIME::HELIUM_ACCUMULATION;
-                if ((utils::Compare(m_Mass, MASS_DOUBLE_DETONATION_CO) >= 0) && (utils::Compare(m_HeShell, WD_HE_SHELL_MCRIT_DETONATION) >= 0)) {
-                    m_HeShellDetonation = true;
-                }
-            }
-        } 
-        else if (utils::Compare(logMdot, logMassTransferCrit) > 0) {
-            regime = ACCRETION_REGIME::HELIUM_OPT_THICK_WINDS;
-        } 
-        else {
-            regime = ACCRETION_REGIME::HELIUM_STABLE_BURNING;
-            if ((utils::Compare(logMdot, COWD_LOG_MDOT_MIN_OFF_CENTER_IGNITION) > 0) && (utils::Compare(m_Mass, COWD_MASS_MIN_OFF_CENTER_IGNITION) > 0)) {
-                m_OffCenterIgnition = true;
-            }
-        }
-    } 
-    else {
-        // The following coefficients in logMassTransfer limits come from quadratic fits to Nomoto+ 2007 results (table 5) in Mass vs log10 Mdot space, to cover the low-mass end.
-        double m_Mass_2 = m_Mass * m_Mass;
-        double logMassTransferCrit   = WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_0 + WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_1 * m_Mass + WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 * m_Mass_2;
-        double logMassTransferStable = WD_LOG_MT_LIMIT_NOMOTO_STABLE_0   + WD_LOG_MT_LIMIT_NOMOTO_STABLE_1   * m_Mass + WD_LOG_MT_LIMIT_NOMOTO_STABLE_2   * m_Mass_2;
-
-        if (utils::Compare(logMdot, logMassTransferStable) < 0) {
-            regime = ACCRETION_REGIME::HYDROGEN_FLASHES;
-        } 
-        else if (utils::Compare(logMdot, logMassTransferCrit) > 0) {
-            regime = ACCRETION_REGIME::HYDROGEN_OPT_THICK_WINDS;
-        } 
-        else {
-            regime = ACCRETION_REGIME::HYDROGEN_STABLE_BURNING;
-        }
-    }
-
-    return regime;
-}
-
-
 /*
  * Specifies next stage, if the star changes its phase.
  *

src/COWD.h

@@ -43,8 +43,6 @@ class COWD: virtual public BaseStar, public WhiteDwarfs {
                                                                                                                                             p_Metallicity, 
                                                                                                                                             WD_Baryon_Number.at(STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF)); }
 
-    ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, const double p_DonorThermalMassLossRate);                                                    // Get the current accretion regime. Can also change flags related to SN events.
-
 protected:
 
     void Initialise() {
@@ -69,12 +67,7 @@ class COWD: virtual public BaseStar, public WhiteDwarfs {
                                                const double p_Metallicity) const    { return CalculateLuminosityOnPhase_Static(p_Mass, p_Time, p_Metallicity); }
     double          CalculateLuminosityOnPhase() const                              { return CalculateLuminosityOnPhase(m_Mass, m_Age, m_Metallicity); }        // Use class member variables
 
-    DBL_DBL         CalculateMassAcceptanceRate(const double p_DonorMassRate,
-                                                const bool   p_IsHeRich);  
-    DBL_DBL         CalculateMassAcceptanceRate(const double p_DonorMassRate,
-                                                const double p_AccretorMassRate,
-                                                const bool   p_IsHeRich)            { return CalculateMassAcceptanceRate(p_DonorMassRate, p_IsHeRich); }        // Ignore the input accretion rate for WDs
-
+
     STELLAR_TYPE    EvolveToNextPhase();
     bool            IsSupernova() const                                             { return m_HeShellDetonation || IsMassAboveChandrasekhar(); };                                             
     bool            ShouldEvolveOnPhase() const                                     { return m_OffCenterIgnition ? false : !IsSupernova(); };                   // From https://ui.adsabs.harvard.edu/abs/2017MNRAS.472.1593W/abstract around the end of section 3.2. Also, allows SN.                                               

src/HG.h

@@ -81,7 +81,7 @@ class HG: virtual public BaseStar, public GiantBranch {
     double          CalculateCoreMassOnPhaseIgnoringPreviousCoreMass(const double p_Mass, const double p_Time) const;                                                           //  Ignore previous core mass constraint when computing expected core mass
 
     double          CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const;
-    double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const        { return 0.25; }                                                                            // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
+    double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const        { return HURLEY_HJELLMING_WEBBINK_QCRIT_HG; }
 
     double          CalculateHeCoreMassAtPhaseEnd() const                           { return m_CoreMass; }                                                                      // McHe(HG) = Core Mass
     double          CalculateHeCoreMassOnPhase() const                              { return m_CoreMass; }                                                                      // McHe(HG) = Core Mass

src/HeGB.h

@@ -59,7 +59,7 @@ class HeGB: virtual public BaseStar, public HeHG {
 
     // member functions - alphabetically
     double      CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const ;
-    double      CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                            { return 1.28; }                                        // From BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
+    double      CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                            { return HURLEY_HJELLMING_WEBBINK_QCRIT_HE_GIANT; }
 
     double      CalculateLuminosityOnPhase(const double p_CoreMass, const double p_GBPB, const double p_GBPD) const { return CalculateLuminosityOnPhase_Static(p_CoreMass, p_GBPB, p_GBPD); }
     double      CalculateLuminosityOnPhase() const                                                                  { return CalculateLuminosityOnPhase(m_CoreMass, m_GBParams[static_cast<int>(GBP::B)], m_GBParams[static_cast<int>(GBP::D)]); }

src/HeHG.h

@@ -69,7 +69,7 @@ class HeHG: virtual public BaseStar, public HeMS {
             double          CalculateCoreMassOnPhase() const                                                        { return m_COCoreMass; }                                                // Mc(HeMS) = McCOMass
 
             double          CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const;
-            double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                { return 1.28; }                                                        // From BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
+            double          CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                { return HURLEY_HJELLMING_WEBBINK_QCRIT_HE_GIANT; }
 
             void            CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams);
             void            CalculateGBParams()                                                                     { CalculateGBParams(m_Mass0, m_GBParams); }                             // Use class member variables