diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 06fc0096a..684da3662 100644 --- a/online-docs/pages/whats-new.rst +++ b/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| diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 92440647a..317ac2f55 100644 --- a/src/BaseBinaryStar.cpp +++ b/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; diff --git a/src/BaseStar.h b/src/BaseStar.h index a340f073d..2852523dc 100644 --- a/src/BaseStar.h +++ b/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 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/COWD.cpp b/src/COWD.cpp index 31b616c42..b923f1cb8 100755 --- a/src/COWD.cpp +++ b/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. * diff --git a/src/COWD.h b/src/COWD.h index 4414eaf93..f3aebab0f 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 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. 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/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/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/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.h b/src/ONeWD.h index a4b757cc6..a0103b5d8 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: @@ -71,6 +70,7 @@ class ONeWD: virtual public BaseStar, public WhiteDwarfs { double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Mass, m_Age, m_Metallicity); } // Use class member variables + STELLAR_TYPE EvolveToNextPhase(); bool IsSupernova() const; bool ShouldEvolveOnPhase() const; diff --git a/src/Remnants.h b/src/Remnants.h index 9bf18d109..2d40afcf1 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/Star.h b/src/Star.h index 60f44fb4f..be87e0d1c 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 494647678..feef739c3 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -1,26 +1,28 @@ #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 * \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) * - * @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 +54,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; @@ -81,7 +83,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 +103,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; @@ -138,12 +140,46 @@ 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 * - * 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). + * 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. * * double CalculateRadiusOnPhase_Static(const double p_Mass) * @@ -157,10 +193,84 @@ 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 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 = 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.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); +} + + +/* + * 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 std::max(NEUTRON_STAR_RADIUS, 0.0115 * std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds))); + return regime; } diff --git a/src/WhiteDwarfs.h b/src/WhiteDwarfs.h index 266c9ff93..e1fcea1cd 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 m_HeShellDetonation and m_OffCenterIgnition flags. + void ResolveShellChange(const double p_AccretedMass); @@ -54,6 +55,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,10 +74,10 @@ 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); } + 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); } diff --git a/src/changelog.h b/src/changelog.h index ca37c2af3..d6ef6cb3a 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1569,7 +1569,17 @@ // - Replaced the name of the KAPIL2024 tides presctiption with KAPIL2025 to match the reference. // - Updated equilibrium and dynamical tides to be consistent with paper. Most notably, corrected all tidal terms to have l=2, and updated the other indices to n and m. // - Added variables for circularization timescale, synchronization timescales (for both stars), and all the tidal ImKnm potential Love numbers to BSE output. +// 03.20.01 SS/NRS - May 26, 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) +// - Moved white dwarf related constants to constants.h (resolves issue #1351) +// - Set merger on unstable RLOF from WD +// + + +const std::string VERSION_STRING = "03.20.01"; -const std::string VERSION_STRING = "03.20.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index a370f24cf..f7da2f0d1 100755 --- a/src/constants.h +++ b/src/constants.h @@ -427,20 +427,55 @@ 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; -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_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. +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). +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