From 51793401a858f3a30b283e29ae71a04ae52ada28 Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Wed, 7 May 2025 18:49:47 +1000 Subject: [PATCH 01/19] Updates to WD accretion to address issue #1384 --- online-docs/pages/whats-new.rst | 5 +++++ src/ONeWD.cpp | 34 +++++++++++++++++++++++++++++++++ src/ONeWD.h | 6 ++++++ src/WhiteDwarfs.cpp | 26 +++++++++++++++---------- src/changelog.h | 8 +++++++- 5 files changed, 68 insertions(+), 11 deletions(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index dd7d0b244..144943e60 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,11 @@ 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.18.05 May 7, 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.18.02 May 1, 2025** * Changed default for Nanjing lambdas to use enhanced lambdas and interpolate in mass and metallicity diff --git a/src/ONeWD.cpp b/src/ONeWD.cpp index 0439a3c9c..f2a4dcdd2 100755 --- a/src/ONeWD.cpp +++ b/src/ONeWD.cpp @@ -1,5 +1,39 @@ #include "ONeWD.h" +/* For ONeWDs, calculate: + * + * (a) the maximum mass acceptance rate of this star, as the accretor, during mass transfer, and + * (b) the retention efficiency parameter + * + * This currently uses the same prescription as for COWDs, but we may consider different + * prescriptions in the future. + * + * 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 (Msun/Myr) + * @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 ONeWD::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); +} + /* * Allow evolution to a new phase (currently, only SN) * diff --git a/src/ONeWD.h b/src/ONeWD.h index a4b757cc6..f10cbb4a1 100755 --- a/src/ONeWD.h +++ b/src/ONeWD.h @@ -71,6 +71,12 @@ class ONeWD: virtual public BaseStar, public WhiteDwarfs { 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; bool ShouldEvolveOnPhase() const; diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 494647678..fb320cacd 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -13,14 +13,14 @@ * * double CalculateEtaH(const double p_MassTransferRate) * - * @param [IN] p_MassTransferRate Mass transfer rate onto the WD surface (Msun/yr) + * @param [IN] p_MassTransferRate Mass transfer rate onto the WD surface (Msun/Myr) * @return Hydrogen accretion efficiency */ double WhiteDwarfs::CalculateEtaH(const double p_MassTransferRate) { double etaH = 0.0; // default return value - double logMassTransferRate = log10(p_MassTransferRate); + double logMassTransferRate = log10(p_MassTransferRate / MYR_TO_YEAR); double m_Mass_2 = m_Mass * m_Mass; // The following coefficients come from quadratic fits to Nomoto+ 2007 results (table 5) in Mass vs log10 Mdot space, to cover the low-mass end. @@ -52,14 +52,14 @@ double WhiteDwarfs::CalculateEtaH(const double p_MassTransferRate) { * * double CalculateEtaHe(const double p_MassTransferRate) * - * @param [IN] p_MassTransferRate Mass transfer rate onto the WD surface (Msun/yr) + * @param [IN] p_MassTransferRate Mass transfer rate onto the WD surface (Msun/Myr) * @return Helium accretion efficiency */ double WhiteDwarfs::CalculateEtaHe(const double p_MassTransferRate) { double etaHe = 1.0; // default return value - so we can have double detonations - double logMassTransferRate = log10(p_MassTransferRate); + double logMassTransferRate = log10(p_MassTransferRate / MYR_TO_YEAR); // The following coefficients in massTransfer limits come from table A1 in Piersanti+ 2014. double logMdotUppHe = WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_0 + WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_1 * m_Mass; @@ -142,8 +142,7 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const /* * Calculate the radius of a white dwarf - good for all types of WD * - * Hurley et al. 2000, eq 91 (from Tout et al. 1997) - * + * Originally from Eggleton 1986, quoted in Verbunt & Rappaport 1988 and Marsh et al. 2004 (eq. 24) * * double CalculateRadiusOnPhase_Static(const double p_Mass) * @@ -157,12 +156,19 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only expected to come up if asking for the core or remnant radius of a giant star - double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); - double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; + const double MP = 5.7E-4; // Constant + const double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); + const double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; + + double MP_Mass = MP / p_Mass; + double MP_Mass_two_thirds = MP_Mass / std::cbrt(MP / p_Mass); - return std::max(NEUTRON_STAR_RADIUS, 0.0115 * std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds))); -} + double First_Factor = std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds)); + double Pre_Second_Factor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass; + double Second_Factor = std::cbrt(Pre_Second_Factor) / Pre_Second_Factor; + return std::max(NEUTRON_STAR_RADIUS, 0.0114 * First_Factor * Second_Factor); +} /* * Increase shell size after mass transfer episode. Hydrogen and helium shells are kept separately. diff --git a/src/changelog.h b/src/changelog.h index c485f0d9a..7190b5631 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1547,7 +1547,13 @@ // - Fix for issue #1380, which appears when the Loveridge binding energy is so high that lambda is rounded off to zero // 03.18.04 IM - May 4, 2025 - Defect repair: // - Added a check to avoid Loveridge lambda becoming zero when the envelope mass is positive but very small +// 03.18.05 SS - May 7, 2025 - Enhancement: +// - Improvements to mass accretion for massive ONe WDs +// - Added ONe::CalculateMassAcceptanceRate +// - Fix units of logMassTransferRate in WhiteDwarfs::CalculateEtaHe and WhiteDwarfs::CalculateEtaH +// - Update white dwarf mass-radius relation (WhiteDwarfs::CalculateRadiusOnPhase_Static) +// -const std::string VERSION_STRING = "03.18.04"; +const std::string VERSION_STRING = "03.18.05"; # endif // __changelog_h__ From 2bc0b28f76346b70eb58da4ef131edb066dafb5c Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia <64709690+nrsegovia@users.noreply.github.com> Date: Thu, 8 May 2025 14:01:56 +1000 Subject: [PATCH 02/19] WD documentation tweaks --- src/WhiteDwarfs.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index fb320cacd..5204a91cf 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -6,10 +6,10 @@ * thresholds logMdotUppH and logMdotLowH. In Claeys+ 2014, the mass transfer rate is * \dot{M}_{tr} and the thresholds are \dot{M}_{cr,H} and \dot{M}_{cr,H}/8, respectively. * - * However, we have used improved thresholds from Nomoto+ 2007, in which the + * However, we have used improved thresholds from Nomoto+ 2007 in which the * lower boundary is \dot{M}_{stable} and the upper boundary is \dot{M}_{RG}. * More precisely, we implemented quadratic fits to the values in Nomoto+ 2007, - * table 5, as described in Rodriguez+ (in prep). + * table 5, as described in the second COMPAS methods paper (in prep). * * double CalculateEtaH(const double p_MassTransferRate) * @@ -142,7 +142,10 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const /* * Calculate the radius of a white dwarf - good for all types of WD * - * Originally from Eggleton 1986, quoted in Verbunt & Rappaport 1988 and Marsh et al. 2004 (eq. 24) + * Originally from Eggleton 1986, quoted in Verbunt & Rappaport 1988 and Marsh et al. 2004 (eq. 24). + * Compared to the Hurley et al. 2000 prescription, the additional factor that includes MP allows + * for the change to a constant density configuration at low masses (e.g., Zapolsky & Salpeter 1969) + * after mass loss episodes. * * double CalculateRadiusOnPhase_Static(const double p_Mass) * From d95df9f83052eaaf74b4573662c0ac8c20478910 Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia <64709690+nrsegovia@users.noreply.github.com> Date: Thu, 8 May 2025 14:18:59 +1000 Subject: [PATCH 03/19] Fixed HURLEY_HJELLMING_WEBBINk qCrit prescription for Remnants --- src/Remnants.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Remnants.h b/src/Remnants.h index 365152636..9a69f05d7 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -42,6 +42,8 @@ class Remnants: virtual public BaseStar, public HeGB { double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta) { return 0.0; } // Should never be called... + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.59; } + void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables From d170eb612fbc52d829cdfcc499a2e67fbed11216 Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia <64709690+nrsegovia@users.noreply.github.com> Date: Thu, 8 May 2025 14:30:01 +1000 Subject: [PATCH 04/19] adapted some variable names to camelCase --- src/WhiteDwarfs.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 5204a91cf..38e24fe0e 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -166,11 +166,11 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { double MP_Mass = MP / p_Mass; double MP_Mass_two_thirds = MP_Mass / std::cbrt(MP / p_Mass); - double First_Factor = std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds)); - double Pre_Second_Factor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass; - double Second_Factor = std::cbrt(Pre_Second_Factor) / Pre_Second_Factor; + double firstFactor = std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds)); + double preSecondFactor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass; + double secondFactor = std::cbrt(preSecondFactor) / preSecondFactor; - return std::max(NEUTRON_STAR_RADIUS, 0.0114 * First_Factor * Second_Factor); + return std::max(NEUTRON_STAR_RADIUS, 0.0114 * firstFactor * secondFactor); } /* From d9e5fafa23a44623dec6f4cc19821b2e98b2feb5 Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Tue, 13 May 2025 17:54:30 +1000 Subject: [PATCH 05/19] HURLEY_HJELLMING_WEBBINK qCrits moved to constants.h and minor cleanup --- src/CHeB.h | 2 +- src/HG.h | 2 +- src/HeGB.h | 2 +- src/HeHG.h | 2 +- src/HeMS.h | 2 +- src/MS_gt_07.h | 2 +- src/MS_lte_07.h | 2 +- src/Remnants.h | 2 +- src/WhiteDwarfs.cpp | 9 ++++----- src/constants.h | 10 +++++++++- 10 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/CHeB.h b/src/CHeB.h index f9a092631..97f486be2 100755 --- a/src/CHeB.h +++ b/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; } diff --git a/src/HG.h b/src/HG.h index 05f48e30e..915872ea3 100755 --- a/src/HG.h +++ b/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 diff --git a/src/HeGB.h b/src/HeGB.h index 93807db3e..56dc03824 100755 --- a/src/HeGB.h +++ b/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(GBP::B)], m_GBParams[static_cast(GBP::D)]); } diff --git a/src/HeHG.h b/src/HeHG.h index 7ba1a210b..9b45a45d5 100755 --- a/src/HeHG.h +++ b/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 diff --git a/src/HeMS.h b/src/HeMS.h index 5ab8939e0..d813084bd 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -89,7 +89,7 @@ class HeMS: virtual public BaseStar, public TPAGB { static double CalculateCoreMass_Luminosity_q_Static(const double p_Mass, const DBL_VECTOR &p_MassCutoffs) { return 3.0; } double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; - 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; } void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams); void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables diff --git a/src/MS_gt_07.h b/src/MS_gt_07.h index 1974137f6..5714d9bbb 100755 --- a/src/MS_gt_07.h +++ b/src/MS_gt_07.h @@ -60,7 +60,7 @@ class MS_gt_07: virtual public BaseStar, public MainSequence { // member functions - alphabetically double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const ; - 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 CalculateMassLossRateHurley(); double CalculateMassTransferRejuvenationFactor(); diff --git a/src/MS_lte_07.h b/src/MS_lte_07.h index b27cb576c..1bbe1ec8f 100755 --- a/src/MS_lte_07.h +++ b/src/MS_lte_07.h @@ -46,7 +46,7 @@ class MS_lte_07: virtual public BaseStar, public MainSequence { // member functions - alphabetically double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const ; - double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.44; } // 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_MS_LTE_07; } double CalculateMassTransferRejuvenationFactor() { return 1.0; } ENVELOPE DetermineEnvelopeType() const { return ENVELOPE::CONVECTIVE; } // Always CONVECTIVE diff --git a/src/Remnants.h b/src/Remnants.h index 9a69f05d7..8d28d674e 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -42,7 +42,7 @@ class Remnants: virtual public BaseStar, public HeGB { double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta) { return 0.0; } // Should never be called... - double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.59; } + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return HURLEY_HJELLMING_WEBBINK_QCRIT_WD; } void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 38e24fe0e..679cec45f 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -143,7 +143,7 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const * Calculate the radius of a white dwarf - good for all types of WD * * Originally from Eggleton 1986, quoted in Verbunt & Rappaport 1988 and Marsh et al. 2004 (eq. 24). - * Compared to the Hurley et al. 2000 prescription, the additional factor that includes MP allows + * Compared to the Hurley et al. 2000 prescription, the additional factor that includes WD_MP allows * for the change to a constant density configuration at low masses (e.g., Zapolsky & Salpeter 1969) * after mass loss episodes. * @@ -159,15 +159,14 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only expected to come up if asking for the core or remnant radius of a giant star - const double MP = 5.7E-4; // Constant const double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); const double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; - double MP_Mass = MP / p_Mass; - double MP_Mass_two_thirds = MP_Mass / std::cbrt(MP / p_Mass); + double MP_Mass = WD_MP / p_Mass; + double MP_Mass_two_thirds = MP_Mass / std::cbrt(WD_MP / p_Mass); double firstFactor = std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds)); - double preSecondFactor = 1 + 3.5 * MP_Mass_two_thirds + MP_Mass; + double preSecondFactor = 1.0 + 3.5 * MP_Mass_two_thirds + MP_Mass; double secondFactor = std::cbrt(preSecondFactor) / preSecondFactor; return std::max(NEUTRON_STAR_RADIUS, 0.0114 * firstFactor * secondFactor); diff --git a/src/constants.h b/src/constants.h index eda1a854d..07faf3fe1 100755 --- a/src/constants.h +++ b/src/constants.h @@ -439,7 +439,15 @@ constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 = -0.98023471; constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.21757267; constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.57319872; constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2137735; - +constexpr double WD_MP = 5.7E-4; + +// Critical mass ratio constants for CalculateCriticalMassRatioHurleyHjellmingWebbink(). +// Based on Hurley+ 2002 section 2.6.1 and BSE code (inverse of the quoted values). +constexpr double HURLEY_HJELLMING_WEBBINK_QCRIT_MS_LTE_07 = 1.44; +constexpr double HURLEY_HJELLMING_WEBBINK_QCRIT_MS_GT_07 = 0.33; +constexpr double HURLEY_HJELLMING_WEBBINK_QCRIT_HG = 0.25; +constexpr double HURLEY_HJELLMING_WEBBINK_QCRIT_HE_GIANT = 1.28; +constexpr double HURLEY_HJELLMING_WEBBINK_QCRIT_WD = 1.59; // coefficients for the calculation of initial angular frequency for Chemically Homogeneous Evolution // Mandel from Butler 2018 From f5b31ffd6fd3088a7f597d2400f0d5504ee3eae8 Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Wed, 14 May 2025 16:15:54 +1000 Subject: [PATCH 06/19] Moved more wd-related constants to constants.h, reorganized DetermineAccretionRegime parameters and properly placed some qCrits in WhiteDwarfs.h --- src/BaseStar.h | 4 ++-- src/COWD.cpp | 8 ++++---- src/COWD.h | 2 +- src/HeWD.cpp | 14 +++++++------- src/HeWD.h | 2 +- src/ONeWD.cpp | 2 +- src/Remnants.h | 2 -- src/Star.h | 4 ++-- src/WhiteDwarfs.cpp | 12 ++++++------ src/WhiteDwarfs.h | 4 +++- src/constants.h | 27 +++++++++++++++++++++++++++ 11 files changed, 54 insertions(+), 27 deletions(-) diff --git a/src/BaseStar.h b/src/BaseStar.h index 9322c0d00..91d84a64c 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -301,8 +301,8 @@ class BaseStar { void ClearCurrentSNEvent() { m_SupernovaDetails.events.current = SN_EVENT::NONE; } // Clear supernova event/state for current timestep void ClearSupernovaStash() { LOGGING->ClearSSESupernovaStash(); } // Clear contents of SSE supernova stash - 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 diff --git a/src/COWD.cpp b/src/COWD.cpp index 31b616c42..9d844caec 100755 --- a/src/COWD.cpp +++ b/src/COWD.cpp @@ -19,7 +19,7 @@ */ DBL_DBL COWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { - m_AccretionRegime = DetermineAccretionRegime(p_IsHeRich, p_DonorMassRate); + m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich); double acceptanceRate = 0.0; // acceptance mass rate - default = 0.0 double fractionAccreted = 0.0; // accretion fraction - default = 0.0 @@ -42,13 +42,13 @@ DBL_DBL COWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bo * * 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) + * ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorMassLossRate, const bool p_HeRich) * - * @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 + * @param [IN] p_HeRich Whether the accreted material is helium-rich or not * @return Current WD accretion regime */ -ACCRETION_REGIME COWD::DetermineAccretionRegime(const bool p_HeRich, const double p_DonorMassLossRate) { +ACCRETION_REGIME COWD::DetermineAccretionRegime(const double p_DonorMassLossRate, const bool p_HeRich) { double logMdot = log10(p_DonorMassLossRate / MYR_TO_YEAR); // logarithm of the accreted mass (M_sun/yr) ACCRETION_REGIME regime = ACCRETION_REGIME::ZERO; diff --git a/src/COWD.h b/src/COWD.h index 4414eaf93..997b59fb7 100755 --- a/src/COWD.h +++ b/src/COWD.h @@ -43,7 +43,7 @@ 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. + ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, const bool p_HeRich); // Get the current accretion regime. Can also change flags related to SN events. protected: diff --git a/src/HeWD.cpp b/src/HeWD.cpp index ef77e3951..9d8f60423 100755 --- a/src/HeWD.cpp +++ b/src/HeWD.cpp @@ -21,7 +21,7 @@ */ DBL_DBL HeWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { - m_AccretionRegime = DetermineAccretionRegime(p_IsHeRich, p_DonorMassRate); // Check if accretion leads to stage switch for WDs and returns retention efficiency as well. + m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich); // Check if accretion leads to stage switch for WDs and returns retention efficiency as well. double acceptanceRate = 0.0; // acceptance mass rate - default = 0.0 double fractionAccreted = m_AccretionRegime == ACCRETION_REGIME::HELIUM_WHITE_DWARF_HYDROGEN_FLASHES ? 0.0 : 1.0; // accretion fraction - default = 1.0, but flashes restrict accumulation @@ -37,19 +37,19 @@ DBL_DBL HeWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bo * * The accretion regime is one of the options listed in enum ACCRETION_REGIME (constants.h) * - * ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, const double p_DonorMassLossRate) + * ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorMassLossRate, const bool p_HeRich) * - * @param [IN] p_HeRich Whether the accreted material is helium-rich or not * @param [IN] p_DonorMassRate Donor mass loss rate, in units of Msol / Myr + * @param [IN] p_HeRich Whether the accreted material is helium-rich or not * @return Current WD accretion regime */ -ACCRETION_REGIME HeWD::DetermineAccretionRegime(const bool p_HeRich, const double p_DonorMassRate) { +ACCRETION_REGIME HeWD::DetermineAccretionRegime(const double p_DonorMassRate, const bool p_HeRich) { double Mdot = p_DonorMassRate / MYR_TO_YEAR; // Accreted mass rate (M_sun/yr) ACCRETION_REGIME regime; if (p_HeRich) { if (utils::Compare(Mdot, HEWD_HE_MDOT_CRIT) <= 0) { regime = ACCRETION_REGIME::HELIUM_WHITE_DWARF_HELIUM_SUB_CHANDRASEKHAR; // Could lead to Sub-Chandrasekhar SN Ia - double massSubCh = -4.0e8 * Mdot + 1.34; // Minimum mass for Sub-Chandrasekhar Mass detonation. Eq 62, Belczynski+ 2008. + double massSubCh = WD_BELCZYNSKI_SN_CONSTANT - WD_BELCZYNSKI_SN_LINEAR * Mdot; // Minimum mass for Sub-Chandrasekhar Mass detonation. Eq 62, Belczynski+ 2008. if (utils::Compare(m_Mass, massSubCh) >= 0 ) { m_IsSubChandrasekharTypeIa = true; } @@ -57,8 +57,8 @@ ACCRETION_REGIME HeWD::DetermineAccretionRegime(const bool p_HeRich, const doubl else { regime = ACCRETION_REGIME::HELIUM_WHITE_DWARF_HELIUM_IGNITION; // Could lift degeneracy and evolve into He MS. Requires minimum mass ! on top of the shell size if (utils::Compare(m_Mass, HEWD_MINIMUM_MASS_IGNITION) >= 0) { - if (utils::Compare(Mdot, 1.64e-6) < 0) { // Accretion limit from eq 61, Belczynski+ 2008. - double mCritHeShell = -7.8e-4 * Mdot + 0.13; // Minimum shell mass of He for detonation. Eq 61, Belczynski+ 2008. This helium should not be burnt, but not implemented this yet. Ruiter+ 2014. + if (utils::Compare(Mdot, WD_BELCZYNSKI_IMMEDIATE_FLASH) < 0) { // Accretion limit from eq 61, Belczynski+ 2008. + double mCritHeShell = WD_BELCZYNSKI_MINIMUM_HE_CONSTANT - WD_BELCZYNSKI_MINIMUM_HE_LINEAR * Mdot;// Minimum shell mass of He for ignition. Eq 61, Belczynski+ 2008. This helium should not be burnt, but not implemented this yet. Ruiter+ 2014. if (utils::Compare(m_HeShell, mCritHeShell) >= 0) { m_ShouldRejuvenate = true; } diff --git a/src/HeWD.h b/src/HeWD.h index 8d990e413..9c73909d2 100755 --- a/src/HeWD.h +++ b/src/HeWD.h @@ -41,7 +41,7 @@ class HeWD: virtual public BaseStar, public WhiteDwarfs { return WhiteDwarfs::CalculateLuminosityOnPhase_Static(p_Mass, p_Time, p_Metallicity, WD_Baryon_Number.at(STELLAR_TYPE::HELIUM_WHITE_DWARF)); } - ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, const double p_DonorThermalMassLossRate); + ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, const bool p_HeRich); protected: diff --git a/src/ONeWD.cpp b/src/ONeWD.cpp index f2a4dcdd2..8f1bc152e 100755 --- a/src/ONeWD.cpp +++ b/src/ONeWD.cpp @@ -21,7 +21,7 @@ */ DBL_DBL ONeWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { - m_AccretionRegime = DetermineAccretionRegime(p_IsHeRich, p_DonorMassRate); + m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich); double acceptanceRate = 0.0; // acceptance mass rate - default = 0.0 double fractionAccreted = 0.0; // accretion fraction - default = 0.0 diff --git a/src/Remnants.h b/src/Remnants.h index 8d28d674e..365152636 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -42,8 +42,6 @@ class Remnants: virtual public BaseStar, public HeGB { double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta) { return 0.0; } // Should never be called... - double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return HURLEY_HJELLMING_WEBBINK_QCRIT_WD; } - void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables diff --git a/src/Star.h b/src/Star.h index c9b6e23a4..aea170f2a 100755 --- a/src/Star.h +++ b/src/Star.h @@ -215,8 +215,8 @@ class Star { void ClearCurrentSNEvent() { m_Star->ClearCurrentSNEvent(); } - ACCRETION_REGIME DetermineAccretionRegime(const bool p_HeRich, - const double p_DonorThermalMassLossRate) { return m_Star->DetermineAccretionRegime(p_HeRich, p_DonorThermalMassLossRate); } + ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, + const bool p_HeRich) { return m_Star->DetermineAccretionRegime(p_DonorThermalMassLossRate, p_HeRich); } ENVELOPE DetermineEnvelopeType() const { return m_Star->DetermineEnvelopeType(); } diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 679cec45f..e516b4de7 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -81,7 +81,7 @@ double WhiteDwarfs::CalculateEtaHe(const double p_MassTransferRate) { /* Calculate accretion efficiency as indicated in Piersanti+ 2014, section A3. Their recipe works - * for specific mass and Mdot values, so a better implementation requires interpolation and + * for specific mass and Mdot values, so a better implementation would require interpolation and * extrapolation (specially towards the low-mass end). Right now, we just adopt a * piece-wise approach. Note that the authors also specify that this is based on the first * strong flash only, but we use it for all episodes. @@ -101,19 +101,19 @@ double WhiteDwarfs::CalculateEtaPTY(const double p_MassTransferRate) { // Limits on each conditional statement come from masses from each model in Piersanti+ 2014. The final etaPTY value is based on table A3. if (utils::Compare(m_Mass, 0.6) <= 0) { - etaPTY = 6.0e-3 + 5.1e-2 * massRate + 8.3e-3 * massRate_2 - 3.317e-4 * massRate_3; + etaPTY = WD_PIERSANTI_M060_G0 + WD_PIERSANTI_M060_G1 * massRate + WD_PIERSANTI_M060_G2 * massRate_2 - WD_PIERSANTI_M060_G3 * massRate_3; } else if (utils::Compare(m_Mass, 0.7) <= 0) { - etaPTY = -3.5e-2 + 7.5e-2 * massRate - 1.8e-3 * massRate_2 + 3.266e-5 * massRate_3; + etaPTY = -WD_PIERSANTI_M070_G0 + WD_PIERSANTI_M070_G1 * massRate - WD_PIERSANTI_M070_G2 * massRate_2 + WD_PIERSANTI_M070_G3 * massRate_3; } else if (utils::Compare(m_Mass, 0.81) <= 0) { - etaPTY = 9.3e-2 + 1.8e-2 * massRate + 1.6e-3 * massRate_2 - 4.111e-5 * massRate_3; + etaPTY = WD_PIERSANTI_M081_G0 + WD_PIERSANTI_M081_G1 * massRate + WD_PIERSANTI_M081_G2 * massRate_2 - WD_PIERSANTI_M081_G3 * massRate_3; } else if (utils::Compare(m_Mass, 0.92) <= 0) { - etaPTY = -7.59e-2 + 1.54e-2 * massRate + 4.0e-4 * massRate_2 - 5.905e-6 * massRate_3; + etaPTY = -WD_PIERSANTI_M092_G0 + WD_PIERSANTI_M092_G1 * massRate + WD_PIERSANTI_M092_G2 * massRate_2 - WD_PIERSANTI_M092_G3 * massRate_3; } else { - etaPTY = -0.323 + 4.1e-2 * massRate - 7.0e-4 * massRate_2 + 4.733e-6 * massRate_3; + etaPTY = -WD_PIERSANTI_M102_G0 + WD_PIERSANTI_M102_G1 * massRate - WD_PIERSANTI_M102_G2 * massRate_2 + WD_PIERSANTI_M102_G3 * massRate_3; } return etaPTY; diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index 266c9ff93..54364e1ba 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -54,6 +54,8 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { const bool p_DonorIsGiant, const double p_DonorThermalMassLossRate, const double p_MassLostByDonor); + + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return HURLEY_HJELLMING_WEBBINK_QCRIT_WD; } double CalculateCOCoreMassOnPhase() const { return m_COCoreMass; } // NO-OP @@ -71,7 +73,7 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { double CalculateEtaPTY(const double p_MassIntakeRate); - double Calculatel0Ritter() const { return (m_Metallicity > 0.01) ? 1995262.3 : 31622.8; } // Luminosity constant which depends on metallicity in Ritter 1999, eq 10 + double Calculatel0Ritter() const { return (m_Metallicity > 0.01) ? L0_RITTER_HIGH_Z : L0_RITTER_LOW_Z; } virtual DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { return std::make_tuple(0.0, 0.0); } diff --git a/src/constants.h b/src/constants.h index 07faf3fe1..dc3ad14f9 100755 --- a/src/constants.h +++ b/src/constants.h @@ -426,6 +426,8 @@ constexpr double HEWD_HE_MDOT_CRIT = 2.0E-8; constexpr double HEWD_MINIMUM_MASS_IGNITION = 0.35; // Minimum mass for HeMS burning constexpr double MASS_DOUBLE_DETONATION_CO = 0.9; // Minimum mass for detonation which would yield something similar to SN Ia. Ruiter+ 2014. constexpr double Q_HYDROGEN_BURNING = 6.4E18 * MSOL_TO_G / (SECONDS_IN_YEAR * LSOL); // 6.4E18 is the energy yield of H burning in erg/g as given in Nomoto+ 2007 (2007ApJ...663.1269N) +constexpr double L0_RITTER_HIGH_Z = 1995262.3; // Luminosity constant which depends on metallicity in Ritter 1999, eq 10 +constexpr double L0_RITTER_LOW_Z = 31622.8; constexpr double WD_HE_SHELL_MCRIT_DETONATION = 0.05; // Minimum shell mass of He for detonation. Should be composed of helium (so, exclude burnt material), but not implemented yet. Ruiter+ 2014. constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_0 = -6.84; constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_1 = 1.349; @@ -440,6 +442,31 @@ constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.21757267; constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.57319872; constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2137735; constexpr double WD_MP = 5.7E-4; +constexpr double WD_BELCZYNSKI_SN_CONSTANT = 1.34; // Constant from Eq 62, Belczynski+ 2008 +constexpr double WD_BELCZYNSKI_SN_LINEAR = 4.0E8; // Linear term factor from Eq 62, Belczynski+ 2008 +constexpr double WD_BELCZYNSKI_IMMEDIATE_FLASH = 1.64E-6; // Accretion limit from eq 61, Belczynski+ 2008. +constexpr double WD_BELCZYNSKI_MINIMUM_HE_CONSTANT = 0.13; // Constant from Eq 61, Belczynski+ 2008 +constexpr double WD_BELCZYNSKI_MINIMUM_HE_LINEAR = 7.8E-4; // Linear term factor from Eq 61, Belczynski+ 2008 +constexpr double WD_PIERSANTI_M060_G0 = 6.0E-3; // This and the following PIERSANTI constants follow table A3 in Piersanti+2014 +constexpr double WD_PIERSANTI_M060_G1 = 5.1E-2; +constexpr double WD_PIERSANTI_M060_G2 = 8.3E-3; +constexpr double WD_PIERSANTI_M060_G3 = 3.317E-4; +constexpr double WD_PIERSANTI_M070_G0 = 3.5E-2; +constexpr double WD_PIERSANTI_M070_G1 = 7.5E-2; +constexpr double WD_PIERSANTI_M070_G2 = 1.8E-3; +constexpr double WD_PIERSANTI_M070_G3 = 3.266E-5; +constexpr double WD_PIERSANTI_M081_G0 = 9.3E-2; +constexpr double WD_PIERSANTI_M081_G1 = 1.8E-2; +constexpr double WD_PIERSANTI_M081_G2 = 1.6E-3; +constexpr double WD_PIERSANTI_M081_G3 = 4.111E-5; +constexpr double WD_PIERSANTI_M092_G0 = 7.59E-2; +constexpr double WD_PIERSANTI_M092_G1 = 1.54E-2; +constexpr double WD_PIERSANTI_M092_G2 = 4.0E-4; +constexpr double WD_PIERSANTI_M092_G3 = 5.905E-6; +constexpr double WD_PIERSANTI_M102_G0 = 3.23E-1; +constexpr double WD_PIERSANTI_M102_G1 = 4.1E-2; +constexpr double WD_PIERSANTI_M102_G2 = 7.0E-4; +constexpr double WD_PIERSANTI_M102_G3 = 4.733E-6; // Critical mass ratio constants for CalculateCriticalMassRatioHurleyHjellmingWebbink(). // Based on Hurley+ 2002 section 2.6.1 and BSE code (inverse of the quoted values). From b87bee8f352729076e9bea51e0f03cfe5a1ef6ef Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Thu, 15 May 2025 17:11:55 +1000 Subject: [PATCH 07/19] Moved DetermineAccretionRegime to WhiteDwarfs --- src/COWD.cpp | 67 ------------------------------------------- src/COWD.h | 2 -- src/ONeWD.h | 1 - src/WhiteDwarfs.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++++ src/WhiteDwarfs.h | 3 +- 5 files changed, 72 insertions(+), 71 deletions(-) diff --git a/src/COWD.cpp b/src/COWD.cpp index 9d844caec..5b4db9803 100755 --- a/src/COWD.cpp +++ b/src/COWD.cpp @@ -33,73 +33,6 @@ DBL_DBL COWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bo } -/* - * 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 double p_DonorMassLossRate, const bool p_HeRich) - * - * @param [IN] p_DonorMassLossRate Donor mass loss rate, in units of Msol / Myr - * @param [IN] p_HeRich Whether the accreted material is helium-rich or not - * @return Current WD accretion regime - */ -ACCRETION_REGIME COWD::DetermineAccretionRegime(const double p_DonorMassLossRate, const bool p_HeRich) { - - 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. * diff --git a/src/COWD.h b/src/COWD.h index 997b59fb7..a90d42f5b 100755 --- a/src/COWD.h +++ b/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 double p_DonorThermalMassLossRate, const bool p_HeRich); // Get the current accretion regime. Can also change flags related to SN events. - protected: void Initialise() { diff --git a/src/ONeWD.h b/src/ONeWD.h index f10cbb4a1..0b7c58420 100755 --- a/src/ONeWD.h +++ b/src/ONeWD.h @@ -41,7 +41,6 @@ class ONeWD: virtual public BaseStar, public WhiteDwarfs { p_Time, p_Metallicity, WD_Baryon_Number.at(STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF)); } - protected: diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index e516b4de7..f4f5f0966 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -1,6 +1,8 @@ #include "WhiteDwarfs.h" #include "NS.h" + + /* Calculate eta_hydrogen from Claeys+ 2014, appendix B. This parameter depends * on three regimes for the mass transfer rate, which here are distinguished by the * thresholds logMdotUppH and logMdotLowH. In Claeys+ 2014, the mass transfer rate is @@ -172,6 +174,74 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { return std::max(NEUTRON_STAR_RADIUS, 0.0114 * firstFactor * secondFactor); } + +/* + * 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 double p_DonorMassLossRate, const bool p_HeRich) + * + * @param [IN] p_DonorMassLossRate Donor mass loss rate, in units of Msol / Myr + * @param [IN] p_HeRich Whether the accreted material is helium-rich or not + * @return Current WD accretion regime + */ +ACCRETION_REGIME WhiteDwarfs::DetermineAccretionRegime(const double p_DonorMassLossRate, const bool p_HeRich) { + + 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; +} + + /* * Increase shell size after mass transfer episode. Hydrogen and helium shells are kept separately. * Only applies the full mass increase from accretion to one of the shells. Does not account for, e.g, diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index 54364e1ba..b2d09d122 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -31,7 +31,8 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::OTHER; } // Not A, B, C, or NONE - + ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, const bool p_HeRich); // Get the current accretion regime. Can also change flags related to SN events. + void ResolveShellChange(const double p_AccretedMass); From 12d5970860952c326c88a04d66eaa0a36e83730c Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Thu, 15 May 2025 17:30:25 +1000 Subject: [PATCH 08/19] small tweak to what's new update --- online-docs/pages/whats-new.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 144943e60..653844d7f 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,10 +3,11 @@ 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.18.05 May 7, 2025** +**03.18.08 May 7, 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. +* Moved white dwarf related constants to constants.h **03.18.02 May 1, 2025** From 87ba9bc81464158de0bed114c611278d3f169a3a Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Thu, 15 May 2025 17:41:06 +1000 Subject: [PATCH 09/19] another minor edit to What's New --- online-docs/pages/whats-new.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 653844d7f..b16d62cdf 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,7 +3,7 @@ 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.18.08 May 7, 2025** +**03.18.08 May 15, 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. From 67d07274d2412877f0e78acdd7fd552a5830c77e Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Mon, 19 May 2025 14:12:00 +1000 Subject: [PATCH 10/19] Moving CalculateMassAcceptanceRate to WhiteDwarfs for COWDs and ONeWDs --- src/COWD.cpp | 33 --------------------------------- src/COWD.h | 7 +------ src/NS.h | 2 ++ src/ONeWD.cpp | 34 ---------------------------------- src/ONeWD.h | 7 +------ src/WhiteDwarfs.cpp | 32 ++++++++++++++++++++++++++++++++ src/WhiteDwarfs.h | 4 ++-- 7 files changed, 38 insertions(+), 81 deletions(-) diff --git a/src/COWD.cpp b/src/COWD.cpp index 5b4db9803..b923f1cb8 100755 --- a/src/COWD.cpp +++ b/src/COWD.cpp @@ -1,38 +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_DonorMassRate, p_IsHeRich); - - 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); -} - - /* * Specifies next stage, if the star changes its phase. * diff --git a/src/COWD.h b/src/COWD.h index a90d42f5b..f3aebab0f 100755 --- a/src/COWD.h +++ b/src/COWD.h @@ -67,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. diff --git a/src/NS.h b/src/NS.h index 0546ea2c1..a789371e6 100755 --- a/src/NS.h +++ b/src/NS.h @@ -101,6 +101,8 @@ class NS: virtual public BaseStar, public Remnants { double CalculateBirthMagneticField(); double CalculateBirthSpinPeriod(); + + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.0; } static double CalculateLuminosityOnPhase_Static(const double p_Mass, const double p_Time); double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase_Static(m_Mass, m_Age); } // Use class member variables diff --git a/src/ONeWD.cpp b/src/ONeWD.cpp index 8f1bc152e..0439a3c9c 100755 --- a/src/ONeWD.cpp +++ b/src/ONeWD.cpp @@ -1,39 +1,5 @@ #include "ONeWD.h" -/* For ONeWDs, calculate: - * - * (a) the maximum mass acceptance rate of this star, as the accretor, during mass transfer, and - * (b) the retention efficiency parameter - * - * This currently uses the same prescription as for COWDs, but we may consider different - * prescriptions in the future. - * - * 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 (Msun/Myr) - * @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 ONeWD::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { - - m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich); - - 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); -} - /* * Allow evolution to a new phase (currently, only SN) * diff --git a/src/ONeWD.h b/src/ONeWD.h index 0b7c58420..a0103b5d8 100755 --- a/src/ONeWD.h +++ b/src/ONeWD.h @@ -70,12 +70,7 @@ class ONeWD: virtual public BaseStar, public WhiteDwarfs { 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; bool ShouldEvolveOnPhase() const; diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index f4f5f0966..feef739c3 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -140,6 +140,38 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const return (635.0 * p_Mass * PPOW(p_Metallicity, 0.4)) / PPOW(p_BaryonNumber * (p_Time + 0.1), 1.4); } +/* Calculate: + * + * (a) the maximum mass acceptance rate of this star, as the accretor, during mass transfer, and + * (b) the retention efficiency parameter + * + * Currently used for COWDs and ONeWDs + * + * 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 WhiteDwarfs::CalculateMassAcceptanceRate(const double p_DonorMassRate, const bool p_IsHeRich) { + + m_AccretionRegime = DetermineAccretionRegime(p_DonorMassRate, p_IsHeRich); + + 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); +} /* * Calculate the radius of a white dwarf - good for all types of WD diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index b2d09d122..08028ef2a 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -76,8 +76,8 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { double Calculatel0Ritter() const { return (m_Metallicity > 0.01) ? L0_RITTER_HIGH_Z : L0_RITTER_LOW_Z; } - virtual DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, - const bool p_IsHeRich) { return std::make_tuple(0.0, 0.0); } + 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); } From 8e5edf2e4bf7f23e1c4ed27968fa7df94c8ae8fc Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Mon, 19 May 2025 14:47:44 +1000 Subject: [PATCH 11/19] rounded Nomoto 2007 fit constants to 4 decimals --- src/constants.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/constants.h b/src/constants.h index dc3ad14f9..d87e4f5fc 100755 --- a/src/constants.h +++ b/src/constants.h @@ -435,12 +435,12 @@ constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_0 = -8.115; constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_1 = 2.29; constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_0 = -8.313; constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_1 = 1.018; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_0 = -8.33017155; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_1 = 2.88247131; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 = -0.98023471; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.21757267; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.57319872; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2137735; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_0 = -8.3302; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_1 = 2.8825; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 = -0.9802; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.2176; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.5732; +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2138; constexpr double WD_MP = 5.7E-4; constexpr double WD_BELCZYNSKI_SN_CONSTANT = 1.34; // Constant from Eq 62, Belczynski+ 2008 constexpr double WD_BELCZYNSKI_SN_LINEAR = 4.0E8; // Linear term factor from Eq 62, Belczynski+ 2008 From 535608a7d1629c94ca458f91d55011065d1e643b Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Mon, 19 May 2025 15:00:24 +1000 Subject: [PATCH 12/19] Added CalculateCriticalMassRatioHurleyHjellmingWebbink to Remnants.h --- src/Remnants.h | 2 ++ src/WhiteDwarfs.h | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Remnants.h b/src/Remnants.h index 365152636..023f091d4 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -42,6 +42,8 @@ class Remnants: virtual public BaseStar, public HeGB { double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate, const double p_massTransferEfficiencyBeta) { return 0.0; } // Should never be called... + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.0; } + void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index 08028ef2a..dca6f6f96 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -56,7 +56,7 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { const double p_DonorThermalMassLossRate, const double p_MassLostByDonor); - double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return HURLEY_HJELLMING_WEBBINK_QCRIT_WD; } + double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return HURLEY_HJELLMING_WEBBINK_QCRIT_WD; } double CalculateCOCoreMassOnPhase() const { return m_COCoreMass; } // NO-OP From 168e81d2a21391c587cde3c7ca9a4cb79e20b885 Mon Sep 17 00:00:00 2001 From: jeffriley Date: Tue, 20 May 2025 08:39:46 +1000 Subject: [PATCH 13/19] Update whats-new.rst --- online-docs/pages/whats-new.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index b16d62cdf..07f8b270d 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -7,7 +7,6 @@ Following is a brief list of important updates to the COMPAS code. A complete r * 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. -* Moved white dwarf related constants to constants.h **03.18.02 May 1, 2025** From 7863d35b8ea5ac378d3e50a12737f38827dfa3df Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Thu, 22 May 2025 18:05:33 +1000 Subject: [PATCH 14/19] Implemented merger for unstable MT from WDs when the HURLEY_HJELLMING_WEBBING qCrit prescription is enabled --- src/BaseBinaryStar.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ceb55ca13..5b3f8d5a6 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2086,6 +2086,9 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // NOTE: Critical mass ratio is defined as mAccretor/mDonor double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate(), m_FractionAccreted); isUnstable = utils::Compare((m_Accretor->Mass() / m_Donor->Mass()), qCrit) < 0; + if ((OPTIONS->QCritPrescription() == QCRIT_PRESCRIPTION::HURLEY_HJELLMING_WEBBINK) && (m_Donor->IsOneOf(WHITE_DWARFS))) { + m_Flags.stellarMerger = true; + } } else { // determine stability based on zetas isUnstable = (utils::Compare(m_ZetaStar, m_ZetaLobe) < 0); From 44a139cea74c10dbc3be207805e749705da8bc7b Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Fri, 23 May 2025 14:54:58 +1000 Subject: [PATCH 15/19] Moved white dwarf merger code block so it is triggered under the right conditions with any qCrit prescription --- src/BaseBinaryStar.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 5b3f8d5a6..942aaf28f 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1652,7 +1652,7 @@ 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; } @@ -2086,9 +2086,6 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // NOTE: Critical mass ratio is defined as mAccretor/mDonor double qCrit = m_Donor->CalculateCriticalMassRatio(m_Accretor->IsDegenerate(), m_FractionAccreted); isUnstable = utils::Compare((m_Accretor->Mass() / m_Donor->Mass()), qCrit) < 0; - if ((OPTIONS->QCritPrescription() == QCRIT_PRESCRIPTION::HURLEY_HJELLMING_WEBBINK) && (m_Donor->IsOneOf(WHITE_DWARFS))) { - m_Flags.stellarMerger = true; - } } else { // determine stability based on zetas isUnstable = (utils::Compare(m_ZetaStar, m_ZetaLobe) < 0); From 3f7ca2f3a0d3cc6515cdec15adff8c5d9cc174d8 Mon Sep 17 00:00:00 2001 From: Nicolas Rodriguez-Segovia Date: Fri, 23 May 2025 14:57:35 +1000 Subject: [PATCH 16/19] Update changelog --- src/changelog.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/changelog.h b/src/changelog.h index 654162a59..01c11312a 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1561,6 +1561,7 @@ // - Fix units of logMassTransferRate in WhiteDwarfs::CalculateEtaHe and WhiteDwarfs::CalculateEtaH // - Update white dwarf mass-radius relation (WhiteDwarfs::CalculateRadiusOnPhase_Static) // - Moved white dwarf related constants to constants.h (resolves issue #1351) +// - Set merger on unstable RLOF from WD // const std::string VERSION_STRING = "03.18.08"; From c38f287077f15c923188c0fe93a17f7b905b5bb7 Mon Sep 17 00:00:00 2001 From: jeffriley Date: Sat, 24 May 2025 06:28:56 +1000 Subject: [PATCH 17/19] Update BaseBinaryStar.cpp --- src/BaseBinaryStar.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 942aaf28f..a27a06629 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1641,8 +1641,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 @@ -1658,7 +1658,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } 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; @@ -1670,7 +1670,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 @@ -1699,7 +1699,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; From 4fb151502646e23d56bf1b1ecd6517de18000793 Mon Sep 17 00:00:00 2001 From: Simon Stevenson Date: Mon, 26 May 2025 11:11:31 +1000 Subject: [PATCH 18/19] Added references to WD constants --- src/constants.h | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/constants.h b/src/constants.h index 998726be3..f7da2f0d1 100755 --- a/src/constants.h +++ b/src/constants.h @@ -430,19 +430,19 @@ constexpr double Q_HYDROGEN_BURNING = 6.4E18 * MSOL_TO_G / ( constexpr double L0_RITTER_HIGH_Z = 1995262.3; // Luminosity constant which depends on metallicity in Ritter 1999, eq 10 constexpr double L0_RITTER_LOW_Z = 31622.8; constexpr double WD_HE_SHELL_MCRIT_DETONATION = 0.05; // Minimum shell mass of He for detonation. Should be composed of helium (so, exclude burnt material), but not implemented yet. Ruiter+ 2014. -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_0 = -6.84; -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_1 = 1.349; -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_0 = -8.115; -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_1 = 2.29; -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_0 = -8.313; -constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_1 = 1.018; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_0 = -8.3302; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_1 = 2.8825; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 = -0.9802; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.2176; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.5732; -constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2138; -constexpr double WD_MP = 5.7E-4; +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_0 = -6.84; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_RG_SS_1 = 1.349; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_0 = -8.115; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SS_MF_1 = 2.29; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_0 = -8.313; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_PIERSANTI_SF_Dt_1 = 1.018; // Constant from Piersanti et al. 2014 see Table A1 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_0 = -8.3302; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_1 = 2.8825; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_REDGIANT_2 = -0.9802; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_0 = -9.2176; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_1 = 3.5732; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_LOG_MT_LIMIT_NOMOTO_STABLE_2 = -1.2138; // Quadratic fit to results from Nomoto et al. 2007 Table 5 +constexpr double WD_MP = 5.7E-4; // White dwarf mass parameter determining where mass-radius relation changes for low-mass white dwarfs constexpr double WD_BELCZYNSKI_SN_CONSTANT = 1.34; // Constant from Eq 62, Belczynski+ 2008 constexpr double WD_BELCZYNSKI_SN_LINEAR = 4.0E8; // Linear term factor from Eq 62, Belczynski+ 2008 constexpr double WD_BELCZYNSKI_IMMEDIATE_FLASH = 1.64E-6; // Accretion limit from eq 61, Belczynski+ 2008. From f8d685648601b0dad7e7204ea2f7be262396ae41 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 26 May 2025 14:05:22 +1000 Subject: [PATCH 19/19] Update comments in WhiteDwarfs.h --- src/WhiteDwarfs.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index dca6f6f96..e1fcea1cd 100644 --- a/src/WhiteDwarfs.h +++ b/src/WhiteDwarfs.h @@ -31,7 +31,7 @@ class WhiteDwarfs: virtual public BaseStar, public Remnants { MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::OTHER; } // Not A, B, C, or NONE - ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, const bool p_HeRich); // Get the current accretion regime. Can also change flags related to SN events. + ACCRETION_REGIME DetermineAccretionRegime(const double p_DonorThermalMassLossRate, const bool p_HeRich); // Get the current accretion regime. Can also change m_HeShellDetonation and m_OffCenterIgnition flags. void ResolveShellChange(const double p_AccretedMass);