diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index 85cc6194e..fd0cc2a56 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -236,6 +236,7 @@ stringChoices: # --initial-mass-function: 'KROUPA' # Default: 'KROUPA' # Options: ['KROUPA','UNIFORM','POWERLAW','SALPETER'] # --LBV-mass-loss-prescription: 'HURLEY_ADD' # Default: 'HURLEY_ADD' # Options: ['BELCZYNSKI','HURLEY','HURLEY_ADD','ZERO','NONE'] # --luminous-blue-variable-prescription: 'HURLEY_ADD' # Default: 'HURLEY_ADD' # Options: ['BELCZYNSKI','HURLEY','HURLEY_ADD','ZERO','NONE'] +# --main-sequence-core-mass-prescription: 'MANDEL' # Default: 'MANDEL' # Options: ['SHIKAUCHI','MANDEL','ZERO'] # --mass-loss-prescription: 'MERRITT2024' # Default: 'MERRITT2024' # Options: ['MERRITT2024','BELCZYNSKI2010','HURLEY','ZERO','NONE'] # --OB-mass-loss: 'VINK2021' # Default: 'VINK2021' # Options: ['KRTICKA2018','BJORKLUND2022','VINK2021','VINK2001','ZERO','NONE'] # --OB-mass-loss-prescription: 'VINK2021' # Default: 'VINK2021' # Options: ['KRTICKA2018','BJORKLUND2022','VINK2021','VINK2001','ZERO','NONE'] diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index feb0710fc..5210923e7 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -795,6 +795,15 @@ Default = 4.2 :ref:`Back to Top ` +**--main-sequence-core-mass-prescription** |br| +Main sequence core mass prescription. |br| +Options: {ZERO, MANDEL, SHIKAUCHI} |br| +``ZERO`` : No core mass treatment, set to zero |br| +``MANDEL`` : The core following case A mass transfer is set equal to the expected core mass of a newly formed HG star + with mass equal to that of the donor, scaled by the fraction of the donor's MS lifetime at mass transfer |br| +``SHIKAUCHI`` : Core mass according to Shikauchi et al. (2024) |br| +Default = MANDEL |br| + **--mass-change-fraction** |br| Approximate desired fractional change in stellar mass on phase when setting SSE and BSE timesteps (applied before ``--timestep--multiplier``). |br| Recommended value is 0.005. |br| @@ -1161,7 +1170,8 @@ Default = MULLERMANDEL If TRUE, preserve a larger donor core mass following case A mass transfer. |br| The core is set equal to the expected core mass of a newly formed HG star with mass equal to that of the donor, scaled by the fraction of the donor's MS lifetime at mass transfer. |br| -Default = TRUE +Default = TRUE |br| +DEPRECATION NOTICE: this option has been deprecated and will soon be removed. Please use ``--main-sequence-core-mass-prescription MANDEL`` in future. **--revised-energy-formalism-nandez-ivanova** |br| Enable revised energy formalism of Nandez & Ivanova. |br| @@ -1418,7 +1428,7 @@ Go to :ref:`the top of this page ` for the full alphabetical **Stellar evolution and winds** --use-mass-loss, --check-photon-tiring-limit, --cool-wind-mass-loss-multiplier, --luminous-blue-variable-prescription, --LBV-mass-loss-prescription ---luminous-blue-variable-multiplier, --mass-loss-prescription, --overall-wind-mass-loss-multiplier, --wolf-rayet-multiplier, +--luminous-blue-variable-multiplier, --main-sequence-core-mass-prescription, --mass-loss-prescription, --overall-wind-mass-loss-multiplier, --wolf-rayet-multiplier, --expel-convective-envelope-above-luminosity-threshold, --luminosity-to-mass-threshold, --OB-mass-loss, --OB-mass-loss-prescription, --RSG-mass-loss, --RSG-mass-loss-prescription, --VMS-mass-loss, --vms-mass-loss-prescription, --WR-mass-loss, --WR-mass-loss-prescription diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 3912bde9c..596f8e811 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,12 @@ 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.12.00 Jan 16, 2025** +* Added convective core mass prescription for main sequence stars from Shikauchi+ (2024), describing how the core mass evolves under mass loss and mass gain. +* New command line option ``--main-sequence-core-mass-prescription`` with arguments ``SHIKAUCHI`` (new prescription), ``MANDEL`` (replaces the functionality of ``--retain-core-mass-during-caseA-mass-transfer``), and ``ZERO`` (core mass set to zero, no treatment). +* Added new luminosity prescription for main sequence stars from Shikauchi+ (2024). +* Added treatment for rejuvenation of main sequence accretors when the new prescription is used. + **03.10.00 Nov 29, 2024** Added functionality to log stellar mergers in the BSE switchlog file. diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index a451b0658..892b15f1b 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1736,8 +1736,22 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { double finalMass = (1.0 - phi) * (mass1 + mass2); double initialHydrogenFraction = m_Star1->InitialHydrogenAbundance(); - double finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction; + double finalHydrogenMass; + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && + (utils::Compare(m_Star1->MZAMS(), SHIKAUCHI_LOWER_MASS_LIMIT) >= 0) && + (utils::Compare(m_Star2->MZAMS(), SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) { + + double coreMass1 = m_Star1->MainSequenceCoreMass(); + double coreMass2 = m_Star2->MainSequenceCoreMass(); + + double coreHeliumMass1 = m_Star1->HeliumAbundanceCore() * coreMass1; + double coreHeliumMass2 = m_Star2->HeliumAbundanceCore() * coreMass2; + + finalHydrogenMass = finalMass * initialHydrogenFraction - coreHeliumMass1 - coreHeliumMass2; + } + else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction; + m_SemiMajorAxis = std::numeric_limits::infinity(); // set separation to infinity to avoid subsequent fake interactions with a massless companion (RLOF, CE, etc.) m_Star1->UpdateAfterMerger(finalMass, finalHydrogenMass); @@ -2088,22 +2102,25 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { } else { // have required mass loss massDiffDonor = -massDiffDonor; // set mass difference - m_Donor->UpdateMinimumCoreMass(); // reset the minimum core mass following case A + m_Donor->UpdateTotalMassLossRate(massDiffDonor / (p_Dt * MYR_TO_YEAR)); // update mass loss rate for MS donor + m_Donor->UpdateMainSequenceCoreMass(p_Dt, massDiffDonor / (p_Dt * MYR_TO_YEAR)); // update core mass for MS donor } } if (!m_CEDetails.CEEnow) { // CE flagged? // no - m_MassTransferTrackerHistory = m_Donor == m_Star1 ? MT_TRACKING::STABLE_1_TO_2_SURV : MT_TRACKING::STABLE_2_TO_1_SURV; // record what happened - for later printing - double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness - + m_MassTransferTrackerHistory = m_Donor == m_Star1 ? MT_TRACKING::STABLE_1_TO_2_SURV : MT_TRACKING::STABLE_2_TO_1_SURV; // record what happened - for later printing + + double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness double omegaDonor_pre_MT = m_Donor->Omega(); // used if full donor envelope is removed + m_Accretor->UpdateTotalMassLossRate(massGainAccretor / (p_Dt * MYR_TO_YEAR)); // update mass gain rate for MS accretor + m_Accretor->UpdateMainSequenceCoreMass(p_Dt, massGainAccretor / (p_Dt * MYR_TO_YEAR)); // update core mass for MS accretor + m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor aFinal = CalculateMassTransferOrbit(m_Donor->Mass(), massDiffDonor, *m_Accretor, m_FractionAccreted); // calculate new orbit - m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index f29ecfd9e..590cf1727 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -129,12 +129,13 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_Dt = DEFAULT_INITIAL_DOUBLE_VALUE; m_Tau = DEFAULT_INITIAL_DOUBLE_VALUE; m_Age = 0.0; // ensure age = 0.0 at construction (rather than default initial value) + m_MainSequenceCoreMass = DEFAULT_INITIAL_DOUBLE_VALUE; m_Mass = m_MZAMS; m_Mass0 = m_MZAMS; - m_MinimumCoreMass = 0.0; m_Luminosity = m_LZAMS; m_Radius = m_RZAMS; m_Temperature = m_TZAMS; + m_TotalMassLossRate = DEFAULT_INITIAL_DOUBLE_VALUE; m_ComponentVelocity = Vector3d(); m_OmegaCHE = CalculateOmegaCHE(m_MZAMS, m_Metallicity); @@ -2709,7 +2710,9 @@ double BaseStar::CalculateMassLossRate() { mDot = mDot * OPTIONS->OverallWindMassLossMultiplier(); // apply overall wind mass loss multiplier } - + + UpdateTotalMassLossRate(-mDot); // update total mass loss rate + return mDot; } @@ -4571,12 +4574,15 @@ STELLAR_TYPE BaseStar::EvolveOnPhase(const double p_DeltaTime) { STELLAR_TYPE stellarType = m_StellarType; if (ShouldEvolveOnPhase()) { // evolve timestep on phase + + UpdateMainSequenceCoreMass(p_DeltaTime, -m_Mdot); // update core mass, relevant for MS stars + m_Tau = CalculateTauOnPhase(); m_COCoreMass = CalculateCOCoreMassOnPhase(); m_CoreMass = CalculateCoreMassOnPhase(); m_HeCoreMass = CalculateHeCoreMassOnPhase(); - + m_Luminosity = CalculateLuminosityOnPhase(); // Calculate abundances @@ -4627,7 +4633,7 @@ STELLAR_TYPE BaseStar::ResolveEndOfPhase(const bool p_ResolveEnvelopeLoss) { if (p_ResolveEnvelopeLoss) stellarType = ResolveEnvelopeLoss(); // if required, resolve envelope loss if it occurs if (stellarType == m_StellarType) { // staying on phase? - + m_Tau = CalculateTauAtPhaseEnd(); m_COCoreMass = CalculateCOCoreMassAtPhaseEnd(); diff --git a/src/BaseStar.h b/src/BaseStar.h index ad4262c37..5607bdaea 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -137,13 +137,13 @@ class BaseStar { double LogMetallicitySigma() const { return m_Log10Metallicity; } // sigma in Hurley+ 2000 double LogMetallicityXi() const { return m_Log10Metallicity - LOG10_ZSOL; } // xi in Hurley+ 2000 double Luminosity() const { return m_Luminosity; } + double MainSequenceCoreMass() const { return m_MainSequenceCoreMass; } double Mass() const { return m_Mass; } double Mass0() const { return m_Mass0; } double MassPrev() const { return m_MassPrev; } ST_VECTOR MassTransferDonorHistory() const { return m_MassTransferDonorHistory; } double Mdot() const { return m_Mdot; } double Metallicity() const { return m_Metallicity; } - double MinimumCoreMass() const { return m_MinimumCoreMass; } double MZAMS() const { return m_MZAMS; } double Omega() const { return m_AngularMomentum / CalculateMomentOfInertiaAU(); } double OmegaCHE() const { return m_OmegaCHE; } @@ -185,6 +185,7 @@ class BaseStar { double Temperature() const { return m_Temperature; } double Time() const { return m_Time; } double Timescale(TIMESCALE p_Timescale) const { return m_Timescales[static_cast(p_Timescale)]; } + double TotalMassLossRate() const { return m_TotalMassLossRate; } double TZAMS() const { return m_TZAMS; } virtual ACCRETION_REGIME WhiteDwarfAccretionRegime() const { return ACCRETION_REGIME::ZERO; } double XExponent() const { return m_XExponent; } @@ -212,7 +213,7 @@ class BaseStar { // member functions - alphabetically - void ApplyMassTransferRejuvenationFactor() { m_Age *= CalculateMassTransferRejuvenationFactor(); } // Apply age rejuvenation factor + void ApplyMassTransferRejuvenationFactor() { m_Age *= CalculateMassTransferRejuvenationFactor(); } // Apply age rejuvenation factor void CalculateBindingEnergies(const double p_CoreMass, const double p_EnvMass, const double p_Radius); @@ -282,7 +283,7 @@ class BaseStar { virtual double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return 0.0; } // Only defined for MS stars - virtual double CalculateRemnantRadius() const { return Radius(); } // relevant for MS stars, over-written for GB stars + virtual double CalculateRemnantRadius() const { return Radius(); } // Relevant for MS stars, over-written for GB stars void CalculateSNAnomalies(const double p_Eccentricity); @@ -343,7 +344,7 @@ class BaseStar { void StashSupernovaDetails(const STELLAR_TYPE p_StellarType, const SSE_SN_RECORD_TYPE p_RecordType = SSE_SN_RECORD_TYPE::DEFAULT) { LOGGING->StashSSESupernovaDetails(this, p_StellarType, p_RecordType); } - virtual double TAMSCoreMass() const { return 0.0; } // except MS stars + virtual double TAMSCoreMass() const { return 0.0; } // Except MS stars virtual void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Default is NO-OP virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP @@ -362,8 +363,9 @@ class BaseStar { const double p_MassGainPerTimeStep, const double p_Epsilon) { } // Default is NO-OP - virtual void UpdateMinimumCoreMass() { } // Only set minimal core mass following Main Sequence mass transfer to MS age fraction of TAMS core mass; default is NO-OP - + virtual void UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { } // Set core mass for Main Sequence stars; default is NO-OP + + virtual void UpdateTotalMassLossRate(const double p_MassLossRate) { m_TotalMassLossRate = p_MassLossRate; } // m_TotalMassLossRate = -m_Mdot in SSE, during a mass transfer episode m_TotalMassLossRate = m_MassLossRateInRLOF // printing functions bool PrintDetailedOutput(const int p_Id, const SSE_DETAILED_RECORD_TYPE p_RecordType) const { @@ -379,7 +381,7 @@ class BaseStar { } bool PrintSwitchLog() const { - return OPTIONS->SwitchLog() ? (LOGGING->ObjectSwitchingPersistence() == OBJECT_PERSISTENCE::PERMANENT ? LOGGING->LogSSESwitchLog(this) : true) : true; // Write record to SSE Switchlog log file + return OPTIONS->SwitchLog() ? (LOGGING->ObjectSwitchingPersistence() == OBJECT_PERSISTENCE::PERMANENT ? LOGGING->LogSSESwitchLog(this) : true) : true; // Write record to SSE Switchlog log file } bool PrintSystemParameters(const SSE_SYSPARMS_RECORD_TYPE p_RecordType = SSE_SYSPARMS_RECORD_TYPE::DEFAULT) const { @@ -434,9 +436,9 @@ class BaseStar { double m_HydrogenAbundanceSurface; // Hydrogen abundance at the surface bool m_LBVphaseFlag; // Flag to know if the star satisfied the conditions, at any point in its evolution, to be considered a Luminous Blue Variable (LBV) double m_Luminosity; // Current luminosity (Lsol) + double m_MainSequenceCoreMass; // Core mass of main sequence stars (Msol) double m_Mass; // Current mass (Msol) double m_Mass0; // Current effective initial mass (Msol) - double m_MinimumCoreMass; // Minimum core mass at end of main sequence (MS stars have no core in the Hurley prescription) double m_MinimumLuminosityOnPhase; // Only required for CHeB stars, but only needs to be calculated once per star double m_Mdot; // Current mass loss rate in winds (Msol per yr) MASS_LOSS_TYPE m_DominantMassLossRate; // Current dominant type of wind mass loss @@ -446,6 +448,7 @@ class BaseStar { double m_Tau; // Relative time double m_Temperature; // Current temperature (Tsol) double m_Time; // Current physical time the star has been evolved (Myr) + double m_TotalMassLossRate; // Current mass loss/gain rate from mass transfer or winds (Msol per yr) // Previous timestep variables double m_DtPrev; // Previous timestep diff --git a/src/CH.h b/src/CH.h index d0659b0dd..77741f7f7 100755 --- a/src/CH.h +++ b/src/CH.h @@ -88,7 +88,7 @@ class CH: virtual public BaseStar, public MS_gt_07 { bool ShouldEvolveOnPhase() const { return m_Age < m_Timescales[static_cast(TIMESCALE::tMS)] && (OPTIONS->OptimisticCHE() || Omega() >= m_OmegaCHE); } // Evolve on CHE phase if age in MS timescale and spinning at least as fast as CHE threshold - void UpdateAgeAfterMassLoss(); + void UpdateAgeAfterMassLoss(); }; diff --git a/src/GiantBranch.h b/src/GiantBranch.h index d07ec4fc5..1bbf5cf28 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -135,7 +135,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence { void UpdateInitialMass() { } // NO-OP for most stellar types - void UpdateMinimumCoreMass() { } // NO-OP for most stellar types + void UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { } // NO-OP for most stellar types }; diff --git a/src/HG.cpp b/src/HG.cpp index 5d028a53e..892377999 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -817,10 +817,15 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const #define massCutoffs(x) m_MassCutoffs[static_cast(MASS_CUTOFF::x)] // for convenience and readability - undefined at end of function #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function - double RTMS = MainSequence::CalculateRadiusAtPhaseEnd(p_Mass, p_RZAMS); - double RGB = GiantBranch::CalculateRadiusOnPhase_Static(p_Mass, m_Luminosity, b); + double RTMS; + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) + RTMS = MainSequence::CalculateRadiusAtPhaseEnd(m_Mass, p_RZAMS); // ensures continuity of stellar tracks when SHIKAUCHI core mass prescription is used + else + RTMS = MainSequence::CalculateRadiusAtPhaseEnd(p_Mass, p_RZAMS); - double rx = RGB; // Hurley sse terminlogy (rx) + double RGB = GiantBranch::CalculateRadiusOnPhase_Static(p_Mass, m_Luminosity, b); + + double rx = RGB; // Hurley sse terminlogy (rx) if (utils::Compare(p_Mass, massCutoffs(MFGB)) > 0) { // mass above threshold for He ignition? // yes diff --git a/src/HG.h b/src/HG.h index cf23a5fc6..6c0c12241 100755 --- a/src/HG.h +++ b/src/HG.h @@ -53,10 +53,10 @@ class HG: virtual public BaseStar, public GiantBranch { m_Age = m_Timescales[static_cast(TIMESCALE::tMS)]; // Set age appropriately // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass - // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) - if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) { - double desiredCoreMass = std::min(m_Mass, MinimumCoreMass()); // desired core mass - m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate + // (only relevant if MANDEL or SHIKAUCHI main sequence core mass prescription is used) + if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MainSequenceCoreMass())) < 0) { + double desiredCoreMass = std::min(m_Mass, MainSequenceCoreMass()); // desired core mass + m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate if (m_Mass0 <= 0.0) { // no root found - no solution for estimated core mass m_Mass0 = m_Mass; // if no root found we keep m_Mass0 equal to the total mass } diff --git a/src/MS_gt_07.h b/src/MS_gt_07.h index 9619eb288..62e685132 100755 --- a/src/MS_gt_07.h +++ b/src/MS_gt_07.h @@ -41,6 +41,15 @@ class MS_gt_07: virtual public BaseStar, public MainSequence { void Initialise() { CalculateTimescales(); // Initialise timescales // Age for MS_GT_07 is carried over from CH stars switching to MS after spinning down, so not set to 0.0 here + + // Initialise core mass, luminosity, radius, and temperature if Shikauchi core mass prescription is used + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) { + m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS); + m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; + m_Luminosity = MainSequence::CalculateLuminosityShikauchi(m_MainSequenceCoreMass, m_InitialHeliumAbundance, m_Age); + m_Radius = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Age, m_RZAMS0); + m_Temperature = BaseStar::CalculateTemperatureOnPhase_Static(m_Luminosity, m_Radius); + } } diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 04b7657fe..4f73b506e 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -25,12 +25,14 @@ * double CalculateHeliumAbundanceCoreOnPhase(const double p_Tau) * * @param [IN] p_Tau Fraction of main sequence lifetime - * * @return Helium abundance in the core (Y_c) */ double MainSequence::CalculateHeliumAbundanceCoreOnPhase(const double p_Tau) const { - double heliumAbundanceCoreMax = 1.0 - m_Metallicity; - return ((heliumAbundanceCoreMax - m_InitialHeliumAbundance) * p_Tau) + m_InitialHeliumAbundance; + + // If SHIKAUCHI core mass prescription is used, core helium abundance is calculated with the core mass + return ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) + ? m_HeliumAbundanceCore + : ((1.0 - m_Metallicity - m_InitialHeliumAbundance) * p_Tau) + m_InitialHeliumAbundance; } @@ -45,11 +47,14 @@ double MainSequence::CalculateHeliumAbundanceCoreOnPhase(const double p_Tau) con * double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) * * @param [IN] p_Tau Fraction of main sequence lifetime - * * @return Hydrogen abundance in the core (X_c) */ double MainSequence::CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const { - return m_InitialHydrogenAbundance * (1.0 - p_Tau); + + // If SHIKAUCHI core mass prescription is used, core helium abundance is calculated with the core mass + return ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) + ? 1.0 - m_HeliumAbundanceCore - m_Metallicity + : m_InitialHydrogenAbundance * (1.0 - p_Tau); } @@ -276,13 +281,18 @@ double MainSequence::CalculateLuminosityAtPhaseEnd(const double p_Mass) const { * * @param [IN] p_Time Time (after ZAMS) in Myr * @param [IN] p_Mass Mass in Msol - * @param [IN] p_LZAMS0 Zero Age Main Sequence (ZAMS) Luminosity + * @param [IN] p_LZAMS Zero Age Main Sequence (ZAMS) Luminosity * @return Luminosity on the Main Sequence as a function of time */ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const double p_Mass, const double p_LZAMS) const { #define a m_AnCoefficients // for convenience and readability - undefined at end of function #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function - + + // If SHIKAUCHI core prescription is used, return Shikauchi luminosity + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) { + return CalculateLuminosityShikauchi(m_MainSequenceCoreMass, m_HeliumAbundanceCore, p_Time); + } + const double epsilon = 0.01; double LTMS = CalculateLuminosityAtPhaseEnd(p_Mass); @@ -310,6 +320,44 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl } +/* + * Calculate luminosity on the Main Sequence when Shikauchi et al. (2024) core mass prescription is used + * + * During core hydrogen burning uses eq (A5) from Shikauchi et al. (2024) + * + * When the Main Sequence hook starts (at age 0.99 * tMS) calculates luminosity that smoothly connects the last point + * of core hydrogen burning with the first point of the HG + * + * double CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore, const double p_Age) + * + * @param [IN] p_CoreMass Main sequence core mass in Msol + * @param [IN] p_HeliumAbundanceCore Central helium fraction + * @param [IN] p_Age Current age in Myr + * @return Luminosity on the Main Sequence as a function of current core mass and central helium fraction + */ +double MainSequence::CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore, const double p_Age) const { + DBL_VECTOR L_COEFFICIENTS = std::get<2>(SHIKAUCHI_COEFFICIENTS); + double logMixingCoreMass = std::log10(p_CoreMass); + double tMS = m_Timescales[static_cast(TIMESCALE::tMS)]; + + double logL = L_COEFFICIENTS[0] * logMixingCoreMass + L_COEFFICIENTS[1] * p_HeliumAbundanceCore + L_COEFFICIENTS[2] * logMixingCoreMass * p_HeliumAbundanceCore + L_COEFFICIENTS[3] * logMixingCoreMass * logMixingCoreMass + L_COEFFICIENTS[4] * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[5] * logMixingCoreMass * logMixingCoreMass * logMixingCoreMass + L_COEFFICIENTS[6] * p_HeliumAbundanceCore * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[7] * logMixingCoreMass * logMixingCoreMass * p_HeliumAbundanceCore + L_COEFFICIENTS[8] * logMixingCoreMass * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[9]; + double luminosity = PPOW(10.0, logL); + + if (utils::Compare(p_Age, 0.99 * tMS) >= 0) { // Star in MS hook? + HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); + double luminosityTAMS = clone->Luminosity(); // Get luminosity from clone (with updated Mass0) + delete clone; clone = nullptr; // Return the memory allocated for the clone + + double ageAtHookStart = 0.99 * tMS; + double luminosityAtHookStart = luminosity; // In the hook, core helium abundance fixed at 1-Z and core mass is not changing + + luminosity = (luminosityAtHookStart * (tMS - p_Age) + luminosityTAMS * (p_Age - ageAtHookStart)) / (tMS - ageAtHookStart); // Linear interpolation + + } + return luminosity; +} + + /////////////////////////////////////////////////////////////////////////////////////// // // // RADIUS CALCULATIONS // @@ -477,7 +525,14 @@ double MainSequence::CalculateRadiusAtPhaseEnd(const double p_Mass, const double double MainSequence::CalculateRadiusOnPhase(const double p_Mass, const double p_Time, const double p_RZAMS) const { #define a m_AnCoefficients // for convenience and readability - undefined at end of function #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function - + + // If SHIKAUCHI core prescription is used, return radius that smoothly connects the beginning of MS hook and the beginning of HG + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) { + double tMS = timescales(tMS); + if (utils::Compare(p_Time, 0.99 * tMS) >= 0) // star in MS hook? + return CalculateRadiusTransitionToHG(p_Mass, p_Time, p_RZAMS); + } + const double epsilon = 0.01; double RTMS = CalculateRadiusAtPhaseEnd(p_Mass, p_RZAMS); @@ -539,11 +594,11 @@ double MainSequence::CalculateRadiusOnPhaseTau(const double p_Mass, const double double deltaR = CalculateDeltaR(p_Mass); double gamma = CalculateGamma(p_Mass); - double mu = std::max(0.5, (1.0 - (0.01 * std::max((a[6] / PPOW(p_Mass, a[7])), (a[8] + (a[9] / PPOW(p_Mass, a[10]))))))); // Hurley et al. 2000, eq 7 - double tHook = mu * tBGB; // Hurley et al. 2000, just after eq 5 + double mu = std::max(0.5, (1.0 - (0.01 * std::max((a[6] / PPOW(p_Mass, a[7])), (a[8] + (a[9] / PPOW(p_Mass, a[10]))))))); // Hurley et al. 2000, eq 7 + double tHook = mu * tBGB; // ibid, just after eq 5 double time = tMS * p_Tau; - double tau1 = std::min(1.0, (time / tHook)); // Hurley et al. 2000, eq 14 - double tau2 = std::max(0.0, std::min(1.0, (time - ((1.0 - epsilon) * tHook)) / (epsilon * tHook))); // Hurley et al. 2000, eq 15 + double tau1 = std::min(1.0, (time / tHook)); // ibid, eq 14 + double tau2 = std::max(0.0, std::min(1.0, (time - ((1.0 - epsilon) * tHook)) / (epsilon * tHook))); // ibid, eq 15 // pow() is slow - use multipliaction where it makes sense double tau_3 = p_Tau * p_Tau * p_Tau; @@ -552,18 +607,44 @@ double MainSequence::CalculateRadiusOnPhaseTau(const double p_Mass, const double double tau1_3 = tau1 * tau1 * tau1; double tau2_3 = tau2 * tau2 * tau2; - double logRMS_RZAMS = alphaR * p_Tau; // Hurley et al. 2000, eq 13, part 1 - logRMS_RZAMS += betaR * tau_10; // Hurley et al. 2000, eq 13, part 2 - logRMS_RZAMS += gamma * tau_40; // Hurley et al. 2000, eq 13, part 3 - logRMS_RZAMS += (log10(RTMS / RZAMS) - alphaR - betaR - gamma) * tau_3; // Hurley et al. 2000, eq 13, part 4 - logRMS_RZAMS -= deltaR * (tau1_3 - tau2_3); // Hurley et al. 2000, eq 13, part 5 + double logRMS_RZAMS = alphaR * p_Tau; // ibid, eq 13, part 1 + logRMS_RZAMS += betaR * tau_10; // ibid, eq 13, part 2 + logRMS_RZAMS += gamma * tau_40; // ibid, eq 13, part 3 + logRMS_RZAMS += (log10(RTMS / RZAMS) - alphaR - betaR - gamma) * tau_3; // ibid, eq 13, part 4 + logRMS_RZAMS -= deltaR * (tau1_3 - tau2_3); // ibid, eq 13, part 5 - return RZAMS * PPOW(10.0, logRMS_RZAMS); // rewrite Hurley et al. 2000, eq 13 for R(t) + return RZAMS * PPOW(10.0, logRMS_RZAMS); // rewrite Hurley et al. 2000, eq 13 for R(t) #undef a } +/* + * Calculate radius on the transition from the Main Sequence to the HG when Shikauchi et al. (2024) core mass prescription is used + * + * SHIKAUCHI core mass prescription cannot be used beyond the MS hook (beyond age 0.99 * tMS), and this function smoothly connects + * the radius between the beginning of the hook and the beginning of the HG + * + * + * double CalculateRadiusShikauchiTransitionToHG(const double p_Mass, const double p_Age, double const p_RZAMS) + + * @param [IN] p_Mass Mass in Msol + * @param [IN] p_Age Age in Myr + * @param [IN] p_RZAMS Zero Age Main Sequence (ZAMS) Radius + * @return Radius on the Main Sequence (for age between tHook and tMS) + */ +double MainSequence::CalculateRadiusTransitionToHG(const double p_Mass, const double p_Age, double const p_RZAMS) const { + HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); + double radiusTAMS = clone->Radius(); // Get radius from clone (with updated Mass0) + delete clone; clone = nullptr; // Return the memory allocated for the clone + + double tMS = m_Timescales[static_cast(TIMESCALE::tMS)]; + double tauAtHookStart = 0.99; + double ageAtHookStart = tauAtHookStart * tMS; + double radiusAtHookStart = CalculateRadiusOnPhaseTau(p_Mass, tauAtHookStart); + + return (radiusAtHookStart * (tMS - p_Age) + radiusTAMS * (p_Age - ageAtHookStart)) / (tMS - ageAtHookStart); // Linear interpolation +} /* @@ -578,16 +659,20 @@ double MainSequence::CalculateRadiusOnPhaseTau(const double p_Mass, const double */ double MainSequence::CalculateRadialExtentConvectiveEnvelope() const { double radiusEnvelope0 = m_Radius; + if ( utils::Compare(m_Mass, 1.25) >= 0) radiusEnvelope0 = 0.0; else if (utils::Compare(m_Mass, 0.35) > 0) { - double radiusM035 = CalculateRadiusAtZAMS(0.35); // uses radius of a 0.35 solar mass star at ZAMS rather than at fractional age Tau, but such low-mass stars only grow by a maximum factor of 1.5 [just above Eq. (10) in Hurley, Pols, Tout (2000), so this is a reasonable approximation - radiusEnvelope0 = radiusM035 * std::sqrt((1.25 - m_Mass) / 0.9); + // uses radius of a 0.35 solar mass star at ZAMS rather than at fractional age Tau, + // but such low-mass stars only grow by a maximum factor of 1.5 + // [just above Eq. (10) in Hurley, Pols, Tout (2000)], so this is a reasonable approximation + radiusEnvelope0 = CalculateRadiusAtZAMS(0.35) * std::sqrt((1.25 - m_Mass) / 0.9); } return radiusEnvelope0 * std::sqrt(std::sqrt(1.0 - m_Tau)); } + /* * Calculate the radial extent of the star's convective core (if it has one) * @@ -642,15 +727,16 @@ double MainSequence::CalculateConvectiveCoreMass() const { return (initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass)); } + /* * Calculate the mass of the convective envelope * * Based on section 7.2 (after Eq. 111) of Hurley, Pols, Tout (2000) * * - * double CalculateConvectiveEnvelopeMass() const + * DBL_DBL CalculateConvectiveEnvelopeMass() const * - * @return Mass of convective envelope in Msol + * @return Tuple containing current mass of convective envelope and mass at ZAMS in Msol */ DBL_DBL MainSequence::CalculateConvectiveEnvelopeMass() const { if (utils::Compare(m_Mass, 1.25) > 0) return std::tuple (0.0, 0.0); @@ -664,6 +750,155 @@ DBL_DBL MainSequence::CalculateConvectiveEnvelopeMass() const { } +/* + * Calculate and update the convective core mass and central helium fraction of a main sequence star that loses + * mass either through winds or case A mass transfer according to Shikauchi et al. (2024) + * + * This function also accounts for mass gain by modeling rejuvenation and updates the initial mixing core mass + * and the helium abundance just outside the core + * + * DBL_DBL CalculateMainSequenceCoreMassShikauchi(const double p_Dt, const double p_MassLossRate) + * + * @param [IN] p_Dt Current time step in Myr + * @param [IN] p_MassLossRate Mass loss/gain rate in Msol yr-1 (negative for mass loss, positive for mass gain) + * @return Tuple containing convective core mass and core helium abundance + */ +DBL_DBL MainSequence::CalculateMainSequenceCoreMassShikauchi(const double p_Dt, const double p_MassLossRate) { + + // get Shikauchi coefficients + DBL_VECTOR ALPHA_COEFFICIENTS = std::get<0>(SHIKAUCHI_COEFFICIENTS); + DBL_VECTOR FMIX_COEFFICIENTS = std::get<1>(SHIKAUCHI_COEFFICIENTS); + DBL_VECTOR L_COEFFICIENTS = std::get<2>(SHIKAUCHI_COEFFICIENTS); + + double fmix = FMIX_COEFFICIENTS[0] + FMIX_COEFFICIENTS[1] * std::exp(-m_MZAMS / FMIX_COEFFICIENTS[2]); // Shikauchi et al. (2024), eq (A3) + auto beta = [&](double mass) { return 1.0 - FMIX_COEFFICIENTS[1] * mass / (FMIX_COEFFICIENTS[2] * fmix) * std::exp(-mass / FMIX_COEFFICIENTS[2]); }; // ibid, eq (A4) + auto alpha = [&](double coreMass) { return PPOW(10.0, std::max(-2.0, ALPHA_COEFFICIENTS[1] * coreMass + ALPHA_COEFFICIENTS[2])) + ALPHA_COEFFICIENTS[0]; }; // ibid, eq (A2) + double g = -0.0044 * m_MZAMS + 0.27; // ibid, eq (A7) + auto delta = [&](double centralHeFraction) { return std::min(PPOW(10.0, -(centralHeFraction - m_InitialHeliumAbundance) / (1.0 - m_InitialHeliumAbundance - m_Metallicity) + g), 1.0); }; // ibid, eq (A6) + + // Use boost adaptive ODE solver + controlled_stepper_type controlled_stepper; + state_type x(3); + x[0] = m_HeliumAbundanceCore; + x[1] = std::log(m_MainSequenceCoreMass); + x[2] = m_Mass; + auto ode = [&](const state_type &x, state_type &dxdt, const double) { + dxdt[0] = CalculateLuminosityShikauchi(std::exp(x[1]), x[0], 0.0) / (Q_CNO * std::exp(x[1])); // Shikauchi et al. (2024), eq (13) + dxdt[1] = -alpha(std::exp(x[1])) / (1.0 - alpha(std::exp(x[1])) * x[0]) * dxdt[0] + beta(x[2]) * delta(x[0]) * p_MassLossRate / (YEAR_TO_MYR * x[2]); // ibid, eq (12) + dxdt[2] = p_MassLossRate; // Mass loss/gain rate + }; + integrate_adaptive(controlled_stepper, ode, x, 0.0, p_Dt, p_Dt/100.0); + + double newMixingCoreMass = std::exp(x[1]); // New mixing core mass + double newCentralHeliumFraction = x[0]; // New central helium fraction + double deltaCoreMass = newMixingCoreMass - m_MainSequenceCoreMass; // Difference in core mass + + if (deltaCoreMass > 0.0) { // If the core grows, we need to account for rejuvenation + if (utils::Compare(newMixingCoreMass, m_InitialMainSequenceCoreMass) < 0) { // New core mass less than initial core mass? + // Common factors + double f1 = m_HeliumAbundanceCoreOut - m_InitialHeliumAbundance; + double f2 = m_MainSequenceCoreMass - m_InitialMainSequenceCoreMass; + double f3 = m_MainSequenceCoreMass + deltaCoreMass; + + // Calculate change in helium abundance just outside the core + double deltaYout = f1 / f2 * deltaCoreMass; + + // Calculate the change in core helium abundance, assuming linear profile between Yc and Y0, and that the the accreted gas has helium fraction Y0 + double deltaY = (m_HeliumAbundanceCoreOut - m_HeliumAbundanceCore) / f3 * deltaCoreMass + 0.5 / f3 * f1 / f2 * deltaCoreMass * deltaCoreMass; + newCentralHeliumFraction = m_HeliumAbundanceCore + deltaY; + m_HeliumAbundanceCoreOut += deltaYout; + } + else { // New core mass greater or equal to the initial core mass? + double deltaCoreMass1 = m_InitialMainSequenceCoreMass - m_MainSequenceCoreMass; // Mass accreted up to the initial core mass + double deltaCoreMass2 = deltaCoreMass - deltaCoreMass1; // Remaining accreted mass + newCentralHeliumFraction = (m_MainSequenceCoreMass * m_HeliumAbundanceCore + 0.5 * (m_HeliumAbundanceCoreOut + m_InitialHeliumAbundance) * deltaCoreMass1 + deltaCoreMass2 * m_InitialHeliumAbundance) / (m_MainSequenceCoreMass + deltaCoreMass); + m_HeliumAbundanceCoreOut = m_InitialHeliumAbundance; + m_InitialMainSequenceCoreMass = newMixingCoreMass; + } + } + else + m_HeliumAbundanceCoreOut = newCentralHeliumFraction; // If core did not grow, Y_out = Y_c + + return std::tuple (newMixingCoreMass, std::min(newCentralHeliumFraction, 1.0 - m_Metallicity)); +} + + +/* + * Calculate the initial convective core mass of a main sequence star using Equation (A3) from Shikauchi et al. (2024), + * also used for calculating core mass after MS merger + * + * double CalculateInitialMainSequenceCoreMass(const double p_MZAMS) + * + * @param [IN] p_MZAMS Mass at ZAMS or after merger in Msol + * @return Mass of the convective core at ZAMS or after merger in Msol + */ +double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const { + DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS); + double fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_MZAMS / fmixCoefficients[2]); + return fmix * p_MZAMS; +} + + +/* + * Update the core mass of a main sequence star that loses mass through winds or Case A mass transfer + * When Shikauchi et al. (2024) core prescription is used, also update the core helium abundance and effective age + * + * + * void UpdateMainSequenceCoreMass(const double p_Dt, const double p_MassLossRate) + * + * @param [IN] p_Dt Current timestep in Myr + * @param [IN] p_MassLossRate Mass loss rate either from stellar winds or mass transfer in Msol yr-1 + */ +void MainSequence::UpdateMainSequenceCoreMass(const double p_Dt, const double p_MassLossRate) { + + double mainSequenceCoreMass = m_MainSequenceCoreMass; // default is no change + double heliumAbundanceCore = m_HeliumAbundanceCore; // default is no change + double age = m_Age; // default is no change + + switch (OPTIONS->MainSequenceCoreMassPrescription()) { + case CORE_MASS_PRESCRIPTION::ZERO: + mainSequenceCoreMass = 0.0; + break; + + case CORE_MASS_PRESCRIPTION::MANDEL: + // Calculate the minimum core mass of a main sequence star that loses mass through Case A mass transfer as the core mass of a TAMS star, scaled by the fractional age. + // Only applied to donors as part of binary evolution, not applied to SSE + if ((OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) && (p_MassLossRate < 0.0) && (utils::Compare(p_MassLossRate, -m_Mdot) != 0)) + mainSequenceCoreMass = std::max(m_MainSequenceCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); + break; + + case CORE_MASS_PRESCRIPTION::SHIKAUCHI: + // Set core mass following Shikauchi et al. (2024) + // MZAMS greater than the limit? SHIKAUCHI prescription valid + if (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0) { + // Only proceed with calculations if star is not in MS hook (Yc < 1-Z), time step is not zero, + // and when the mass loss rate argument is equal to the total mass loss rate + // (i.e. total mass loss rate was updated, this prevents the calculation in SSE if it was executed as part of BSE for the same time step) + if ((utils::Compare(m_HeliumAbundanceCore, 1.0 - m_Metallicity) < 0) && (utils::Compare(p_Dt, 0.0) != 0) && (utils::Compare(p_MassLossRate, m_TotalMassLossRate) == 0)) { + + std::tie(mainSequenceCoreMass, heliumAbundanceCore) = CalculateMainSequenceCoreMassShikauchi(p_Dt, p_MassLossRate); // calculate and update the core mass and central helium fraction + + double tMS = m_Timescales[static_cast(TIMESCALE::tMS)]; + age = (heliumAbundanceCore - m_InitialHeliumAbundance) / m_InitialHydrogenAbundance * 0.99 * tMS; // update the effective age based on central helium fraction + } + } + // MZAMS less than the limit? MANDEL prescription used + else { + // Only applied to donors as part of binary evolution, not applied to SSE + if ((p_MassLossRate < 0.0) && (utils::Compare(p_MassLossRate, -m_Mdot) != 0)) + mainSequenceCoreMass = std::max(m_MainSequenceCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); + } + break; + + default: break; // do nothing + } + + m_MainSequenceCoreMass = mainSequenceCoreMass; // update core mass + m_HeliumAbundanceCore = heliumAbundanceCore; // update core helium abundance + m_Age = age; // update age +} + + /////////////////////////////////////////////////////////////////////////////////////// // // // LIFETIME / AGE CALCULATIONS // @@ -733,11 +968,12 @@ void MainSequence::UpdateAgeAfterMassLoss() { double tMS = m_Timescales[static_cast(TIMESCALE::tMS)]; double tBGBprime = CalculateLifetimeToBGB(m_Mass); double tMSprime = MainSequence::CalculateLifetimeOnPhase(m_Mass, tBGBprime); - + m_Age *= tMSprime / tMS; CalculateTimescales(m_Mass, m_Timescales); // must update timescales } + /////////////////////////////////////////////////////////////////////////////////////// // // // MISCELLANEOUS FUNCTIONS / CONTROL FUNCTIONS // @@ -824,23 +1060,6 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_Force) { } -/* - * Update the minimum core mass of a main sequence star that loses mass through Case A mass transfer by - * setting it equal to the core mass of a TAMS star, scaled by the fractional age. - * - * The minimum core mass of the star is updated only if the retain-core-mass-during-caseA-mass-transfer - * option is specified, otherwise it is left unchanged. - * - * - * void UpdateMinimumCoreMass() - * - */ -void MainSequence::UpdateMinimumCoreMass() { - if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { - m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); // update minimum core mass - } -} - /* * Return the expected core mass at terminal age main sequence, i.e., at the start of the HG phase * @@ -867,7 +1086,7 @@ double MainSequence::TAMSCoreMass() const { /* - * Sets the mass and age of a merge product of two main sequence stars + * Sets the mass and age of a merger product of two main sequence stars * (note: treats merger products as main-sequence stars, not CHE, and does not check for MS_lte_07) * * Uses prescription of Wang et al., 2022, https://www.nature.com/articles/s41550-021-01597-5 @@ -881,26 +1100,34 @@ double MainSequence::TAMSCoreMass() const { void MainSequence::UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function - m_Mass = p_Mass; - m_Mass0 = m_Mass; - m_MinimumCoreMass = 0.0; + m_Mass = p_Mass; + m_Mass0 = m_Mass; + m_MainSequenceCoreMass = 0.0; double initialHydrogenFraction = m_InitialHydrogenAbundance; CalculateTimescales(); CalculateGBParams(); - m_Tau = (initialHydrogenFraction - p_HydrogenMass / m_Mass) / initialHydrogenFraction; // assumes uniformly mixed merger product and a uniform rate of H fusion on main sequence + m_Tau = (initialHydrogenFraction - p_HydrogenMass / m_Mass) / initialHydrogenFraction; // assumes uniformly mixed merger product and a uniform rate of H fusion on main sequence m_Age = m_Tau * timescales(tMS); + m_HeliumAbundanceCore = 1.0 - m_Metallicity - p_HydrogenMass / p_Mass; + + m_HydrogenAbundanceCore = 1.0 - m_Metallicity - m_HeliumAbundanceCore; + + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::SHIKAUCHI) && (utils::Compare(m_MZAMS, SHIKAUCHI_LOWER_MASS_LIMIT) >= 0)) { + m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass); // update initial mixing core mass + m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; // update core mass + } + UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); #undef timescales } - /* * Interpolate Ge+ Critical Mass Ratios, for H-rich stars * @@ -1061,3 +1288,71 @@ double MainSequence::InterpolateGeEtAlQCrit(const QCRIT_PRESCRIPTION p_qCritPres return qCritPerMetallicity[1] + (m_Log10Metallicity - logZhi)*(qCritPerMetallicity[1] - qCritPerMetallicity[0])/(logZhi - logZlo); } + + +/* + * Linear interpolation/extrapolation for coefficients from Shikauchi et al. (2024), used for core mass calculations + * + * std::tuple MainSequence::InterpolateShikauchiCoefficients(const double p_Metallicity) + * + * @param [IN] p_Metallicity Metallicity + * @return Tuple containing vectors of coefficients for the specified metallicity + */ +std::tuple MainSequence::InterpolateShikauchiCoefficients(const double p_Metallicity) const { + + DBL_VECTOR alphaCoeff(3, 0.0); + DBL_VECTOR fmixCoeff(3, 0.0); + DBL_VECTOR lCoeff(10, 0.0); + + // Skip calculation if SHIKAUCHI core prescription is not used + if (OPTIONS->MainSequenceCoreMassPrescription() != CORE_MASS_PRESCRIPTION::SHIKAUCHI) + return std::tuple (alphaCoeff, fmixCoeff, lCoeff); + + double logZ = std::log10(p_Metallicity); + + // Coefficients are given for these metallicities + double low = std::log10(0.1 * ZSOL_ASPLUND); + double middle = std::log10(1.0 / 3.0 * ZSOL_ASPLUND); + double high = std::log10(ZSOL_ASPLUND); + + // common factors + double middle_logZ = middle - logZ; + double middle_low = middle - low; + double high_logZ = high - logZ; + double high_middle = high - middle; + double logZ_low = logZ - low; + double logZ_middle = logZ - middle; + + // Linear extrapolation (constant) for metallicity lower than the lowest bound + if (utils::Compare(logZ, low) <= 0) { + alphaCoeff = SHIKAUCHI_ALPHA_COEFFICIENTS[0]; + fmixCoeff = SHIKAUCHI_FMIX_COEFFICIENTS[0]; + lCoeff = SHIKAUCHI_L_COEFFICIENTS[0]; + } + // Linear interpolation between metallicity low and middle + else if ((utils::Compare(logZ, low) > 0) && (utils::Compare(logZ, middle) <= 0)) { + for (size_t i = 0; i < 3; i++) { + alphaCoeff[i] = (SHIKAUCHI_ALPHA_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_ALPHA_COEFFICIENTS[1][i] * logZ_low) / middle_low; + fmixCoeff[i] = (SHIKAUCHI_FMIX_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_FMIX_COEFFICIENTS[1][i] * logZ_low) / middle_low; + } + for (size_t i = 0; i < 10; i++) + lCoeff[i] = (SHIKAUCHI_L_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_L_COEFFICIENTS[1][i] * logZ_low) / middle_low; + } + // Linear interpolation between metallicity middle and high + else if ((utils::Compare(logZ, middle) > 0) && (utils::Compare(logZ, high) < 0)) { + for (size_t i = 0; i < 3; i++) { + alphaCoeff[i] = (SHIKAUCHI_ALPHA_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_ALPHA_COEFFICIENTS[2][i] * logZ_middle) / high_middle; + fmixCoeff[i] = (SHIKAUCHI_FMIX_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_FMIX_COEFFICIENTS[2][i] * logZ_middle) / high_middle; + } + for (size_t i = 0; i < 10; i++) + lCoeff[i] = (SHIKAUCHI_L_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_L_COEFFICIENTS[2][i] * logZ_middle) / high_middle; + } + // Linear extrapolation (constant) for metallicity equal to solar or higher + else { + alphaCoeff = SHIKAUCHI_ALPHA_COEFFICIENTS[2]; + fmixCoeff = SHIKAUCHI_FMIX_COEFFICIENTS[2]; + lCoeff = SHIKAUCHI_L_COEFFICIENTS[2]; + } + + return std::tuple (alphaCoeff, fmixCoeff, lCoeff); +} diff --git a/src/MainSequence.h b/src/MainSequence.h index 36aa11b2c..c035345e5 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -19,10 +19,15 @@ class MainSequence: virtual public BaseStar { MainSequence(const BaseStar& p_BaseStar) : BaseStar(p_BaseStar) {} MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::A; } // Always case A - + + const std::tuple SHIKAUCHI_COEFFICIENTS = InterpolateShikauchiCoefficients(m_Metallicity); // Interpolate Shikauchi coefficients for the given metallicity protected: - + + // member variables + + double m_HeliumAbundanceCoreOut = m_InitialHeliumAbundance; // Helium abundance just outside the core, used for rejuvenation calculations + double m_InitialMainSequenceCoreMass = 0.0; // Initial mass of the mixing core is initialised in MS_gt_07 class // member functions - alphabetically double CalculateAlphaL(const double p_Mass) const; @@ -44,7 +49,7 @@ class MainSequence: virtual public BaseStar { double CalculateCOCoreMassAtPhaseEnd() const { return CalculateCOCoreMassOnPhase(); } // Same as on phase double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(MS) = 0.0 - double CalculateCoreMassAtPhaseEnd() const { return OPTIONS->RetainCoreMassDuringCaseAMassTransfer() ? MinimumCoreMass() : 0.0; } // Accounts for minimal core mass built up prior to mass loss through mass transfer + double CalculateCoreMassAtPhaseEnd() const { return (OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::MANDEL) ? MainSequenceCoreMass() : 0.0; } // Accounts for minimal core mass built up prior to mass loss through mass transfer double CalculateCoreMassOnPhase() const { return 0.0; } // Mc(MS) = 0.0 (Hurley et al. 2000, just before eq 28) double CalculateHeCoreMassAtPhaseEnd() const { return CalculateCoreMassAtPhaseEnd(); } // Same as He core mass @@ -70,7 +75,9 @@ class MainSequence: virtual public BaseStar { double CalculateLuminosityAtPhaseEnd() const { return CalculateLuminosityAtPhaseEnd(m_Mass0); } // Use class member variables double CalculateLuminosityOnPhase(const double p_Time, const double p_Mass, const double p_LZAMS) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0); } // Use class member variables - + double CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore, const double p_Age) const; + DBL_DBL CalculateMainSequenceCoreMassShikauchi(const double p_Dt, const double p_MassLossRate); + double CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const; double CalculateMomentOfInertia() const { return (0.1 * (m_Mass) * m_Radius * m_Radius); } // k2 = 0.1 as defined in Hurley et al. 2000, after eq 109 double CalculatePerturbationMu() const { return 5.0; } // mu(MS) = 5.0 (Hurley et al. 2000, eqs 97 & 98) @@ -83,7 +90,8 @@ class MainSequence: virtual public BaseStar { double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_RZAMS) const; double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass, m_RZAMS); } // Use class member variables double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass, m_Age, m_RZAMS0); } // Use class member variables - + double CalculateRadiusTransitionToHG(const double p_Mass, const double p_Age, double const p_RZAMS) const; + double CalculateTauAtPhaseEnd() const { return 1.0; } // tau = 1.0 at end of MS double CalculateTauOnPhase() const; @@ -98,8 +106,10 @@ class MainSequence: virtual public BaseStar { STELLAR_TYPE EvolveToNextPhase() { return STELLAR_TYPE::HERTZSPRUNG_GAP; } double InterpolateGeEtAlQCrit(const QCRIT_PRESCRIPTION p_qCritPrescription, - const double p_massTransferEfficiencyBeta); // RTW do I need a const here? - + const double p_massTransferEfficiencyBeta); // RTW do I need a const here? + + std::tuple InterpolateShikauchiCoefficients(const double p_Metallicity) const; + bool IsEndOfPhase() const { return !ShouldEvolveOnPhase(); } // Phase ends when age at or after MS timescale void PerturbLuminosityAndRadius() { } // NO-OP @@ -115,8 +125,8 @@ class MainSequence: virtual public BaseStar { void UpdateAfterMerger(double p_Mass, double p_HydrogenMass); void UpdateAgeAfterMassLoss(); // Per Hurley et al. 2000, section 7.1 - - void UpdateMinimumCoreMass(); // Set minimal core mass following Main Sequence mass transfer to MS age fraction of TAMS core mass + + void UpdateMainSequenceCoreMass(const double p_Dt, const double p_MassLossRate); }; diff --git a/src/Options.cpp b/src/Options.cpp index 2b81103d6..c53199418 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -408,6 +408,10 @@ void Options::OptionValues::Initialise() { m_WolfRayetFactor = 1.0; m_ScaleTerminalWindVelocityWithMetallicityPower = 0.0; + // Core mass prescription + m_MainSequenceCoreMassPrescription.type = CORE_MASS_PRESCRIPTION::MANDEL; + m_MainSequenceCoreMassPrescription.typeString = CORE_MASS_PRESCRIPTION_LABEL.at(m_MainSequenceCoreMassPrescription.type); + // Mass transfer options m_UseMassTransfer = true; m_CirculariseBinaryDuringMassTransfer = true; @@ -886,7 +890,7 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt po::value(&p_Options->m_Quiet)->default_value(p_Options->m_Quiet)->implicit_value(true), ("Suppress printing (default = " + std::string(p_Options->m_Quiet ? "TRUE" : "FALSE") + ")").c_str() ) - + ( "retain-core-mass-during-caseA-mass-transfer", po::value(&p_Options->m_RetainCoreMassDuringCaseAMassTransfer)->default_value(p_Options->m_RetainCoreMassDuringCaseAMassTransfer)->implicit_value(true), @@ -1787,9 +1791,14 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt po::value(&p_Options->m_LBVMassLossPrescription.typeString)->default_value(p_Options->m_LBVMassLossPrescription.typeString), ("LBV Mass loss prescription (" + AllowedOptionValuesFormatted("LBV-mass-loss-prescription") + ", default = '" + p_Options->m_LBVMassLossPrescription.typeString + "')").c_str() ) - ( - "mass-loss-prescription", + "main-sequence-core-mass-prescription", + po::value(&p_Options->m_MainSequenceCoreMassPrescription.typeString)->default_value(p_Options->m_MainSequenceCoreMassPrescription.typeString), + + ("Main Sequence core mass prescription (" + AllowedOptionValuesFormatted("main-sequence-core-mass-prescription") + ", default = '" + p_Options->m_MainSequenceCoreMassPrescription.typeString + "')").c_str() + ) + ( + "mass-loss-prescription", po::value(&p_Options->m_MassLossPrescription.typeString)->default_value(p_Options->m_MassLossPrescription.typeString), ("Mass loss prescription (" + AllowedOptionValuesFormatted("mass-loss-prescription") + ", default = '" + p_Options->m_MassLossPrescription.typeString + "')").c_str() ) @@ -2225,7 +2234,12 @@ std::string Options::OptionValues::CheckAndSetOptions() { std::tie(found, m_LBVMassLossPrescription.type) = utils::GetMapKey(m_LBVMassLossPrescription.typeString, LBV_MASS_LOSS_PRESCRIPTION_LABEL, m_LBVMassLossPrescription.type); COMPLAIN_IF(!found, "Unknown LBV Mass Loss Prescription"); } - + + if (!DEFAULTED("main-sequence-core-mass-prescription")) { // main sequence core mass prescription + std::tie(found, m_MainSequenceCoreMassPrescription.type) = utils::GetMapKey(m_MainSequenceCoreMassPrescription.typeString, CORE_MASS_PRESCRIPTION_LABEL, m_MainSequenceCoreMassPrescription.type); + COMPLAIN_IF(!found, "Unknown Main Sequence Core Mass Prescription"); + } + if (!DEFAULTED("mass-loss-prescription")) { // mass loss prescription std::tie(found, m_MassLossPrescription.type) = utils::GetMapKey(m_MassLossPrescription.typeString, MASS_LOSS_PRESCRIPTION_LABEL, m_MassLossPrescription.type); COMPLAIN_IF(!found, "Unknown Mass Loss Prescription"); @@ -2581,6 +2595,7 @@ std::vector Options::AllowedOptionValues(const std::string p_Option case _("kick-magnitude-distribution") : POPULATE_RET(KICK_MAGNITUDE_DISTRIBUTION_LABEL); break; case _("logfile-type") : POPULATE_RET(LOGFILETYPELabel); break; case _("LBV-mass-loss-prescription") : POPULATE_RET(LBV_MASS_LOSS_PRESCRIPTION_LABEL); break; + case _("main-sequence-core-mass-prescription") : POPULATE_RET(CORE_MASS_PRESCRIPTION_LABEL); break; case _("mass-loss-prescription") : POPULATE_RET(MASS_LOSS_PRESCRIPTION_LABEL); break; case _("mass-ratio-distribution") : POPULATE_RET(MASS_RATIO_DISTRIBUTION_LABEL); break; case _("mass-transfer-accretion-efficiency-prescription") : POPULATE_RET(MT_ACCRETION_EFFICIENCY_PRESCRIPTION_LABEL); break; diff --git a/src/Options.h b/src/Options.h index 345bb75e3..495c7ea62 100755 --- a/src/Options.h +++ b/src/Options.h @@ -211,16 +211,17 @@ class Options { // is only shown once per COMPAS run). std::vector> deprecatedOptionStrings = { - { "black-hole-kicks", "black-hole-kicks-mode", false }, - { "chemically-homogeneous-evolution", "chemically-homogeneous-evolution-mode", false }, - { "kick-direction", "kick-direction-distribution", false }, - { "luminous-blue-variable-prescription", "LBV-mass-loss-prescription", false }, - { "mass-transfer", "use-mass-transfer", false }, - { "mass-transfer-thermal-limit-accretor", "mass-transfer-thermal-limit-accretor-multiplier", false }, - { "OB-mass-loss", "OB-mass-loss-prescription", false }, - { "RSG-mass-loss", "RSG-mass-loss-prescription", false }, - { "VMS-mass-loss", "VMS-mass-loss-prescription", false }, - { "WR-mass-loss", "WR-mass-loss-prescription", false } + { "black-hole-kicks", "black-hole-kicks-mode", false }, + { "chemically-homogeneous-evolution", "chemically-homogeneous-evolution-mode", false }, + { "kick-direction", "kick-direction-distribution", false }, + { "luminous-blue-variable-prescription", "LBV-mass-loss-prescription", false }, + { "mass-transfer", "use-mass-transfer", false }, + { "mass-transfer-thermal-limit-accretor", "mass-transfer-thermal-limit-accretor-multiplier", false }, + { "OB-mass-loss", "OB-mass-loss-prescription", false }, + { "retain-core-mass-during-caseA-mass-transfer", "", false }, + { "RSG-mass-loss", "RSG-mass-loss-prescription", false }, + { "VMS-mass-loss", "VMS-mass-loss-prescription", false }, + { "WR-mass-loss", "WR-mass-loss-prescription", false } }; std::vector> deprecatedOptionValues = { @@ -618,6 +619,7 @@ class Options { "logfile-system-parameters-record-types", "logfile-type", + "main-sequence-core-mass-prescription", "mass-change-fraction", "mass-loss-prescription", "mass-ratio-distribution", @@ -991,8 +993,10 @@ class Options { bool m_ExpelConvectiveEnvelopeAboveLuminosityThreshold; // Whether to expel the convective envelope in a pulsation when log_10(L/M) reaches the threshold defined by m_LuminosityToMassThreshold double m_LuminosityToMassThreshold; // Threshold value of log_10(L/M) above which the convective envelope is expelled in a pulsation - + bool m_RetainCoreMassDuringCaseAMassTransfer; // Whether to retain the approximate core mass of a case A donor as a minimum core at end of MS or HeMS (default = false) + + ENUM_OPT m_MainSequenceCoreMassPrescription; // Which MS core prescription ENUM_OPT m_CaseBBStabilityPrescription; // Which prescription for the stability of case BB/BC mass transfer @@ -1502,6 +1506,8 @@ class Options { double LuminousBlueVariableFactor() const { return OPT_VALUE("luminous-blue-variable-multiplier", m_LuminousBlueVariableFactor, true); } LBV_MASS_LOSS_PRESCRIPTION LBVMassLossPrescription() const { return OPT_VALUE("LBV-mass-loss-prescription", m_LBVMassLossPrescription.type, true); } + CORE_MASS_PRESCRIPTION MainSequenceCoreMassPrescription() const { return OPT_VALUE("main-sequence-core-mass-prescription", m_MainSequenceCoreMassPrescription.type, true); } + double MassChangeFraction() const { return m_CmdLine.optionValues.m_MassChangeFraction; } MASS_LOSS_PRESCRIPTION MassLossPrescription() const { return OPT_VALUE("mass-loss-prescription", m_MassLossPrescription.type, true); } diff --git a/src/Star.h b/src/Star.h index c486cc1a9..cfc7a65f9 100755 --- a/src/Star.h +++ b/src/Star.h @@ -118,6 +118,7 @@ class Star { double LambdaKruckow() const { return m_Star->LambdaKruckow(); } double LambdaDewi() const { return m_Star->LambdaDewi(); } double Luminosity() const { return m_Star->Luminosity(); } + double MainSequenceCoreMass() const { return m_Star->MainSequenceCoreMass(); } double Mass() const { return m_Star->Mass(); } double Mass0() const { return m_Star->Mass0(); } double MassPrev() const { return m_Star->MassPrev(); } @@ -142,6 +143,7 @@ class Star { double Tau() const { return m_Star->Tau(); } double Temperature() const { return m_Star->Temperature(); } double Timescale(TIMESCALE p_Timescale) const { return m_Star->Timescale(p_Timescale); } + double TotalMassLossRate() const { return m_Star->TotalMassLossRate(); } double XExponent() const { return m_Star->XExponent(); } @@ -293,10 +295,12 @@ class Star { p_MassGainPerTimeStep, p_Epsilon);} - void UpdateMinimumCoreMass() { m_Star->UpdateMinimumCoreMass(); } - + void UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { m_Star->UpdateMainSequenceCoreMass(p_Dt, p_TotalMassLossRate); } + void UpdatePreviousTimestepDuration() { m_Star->UpdatePreviousTimestepDuration(); } + void UpdateTotalMassLossRate(const double p_MassLossRate) { m_Star->UpdateTotalMassLossRate(p_MassLossRate); } + ACCRETION_REGIME WhiteDwarfAccretionRegime() const { return m_Star->WhiteDwarfAccretionRegime(); } private: diff --git a/src/changelog.h b/src/changelog.h index 255cc8f4d..d977012a1 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1435,7 +1435,13 @@ // - At CHE initialization, stellar spin is set to orbital frequency, unless rotational frequency has been specified by user. This process does not conserve angular momentum (implicitly assuming spin-up in the pre-ZAMS phase). // - When checking for CHE, compare threshold frequency against orbit rather than stellar spin, in case the star has zero frequency (no tides, no user-specified value). // - Moved all CHE rotation related code to ProcessTides(), ensuring that any spin up during binary evolution conserves total angular momentum. +// 03.12.00 AB - Jan 16, 2025 - Enhancement: +// - Added Shikauchi et al. (2024) core mass prescription, describing convective core evolution under mass loss/gain +// - New options: --main-sequence-core-mass-prescription SHIKAUCHI (new prescription), MANDEL (replaces --retain-core-mass-during-caseA-mass-transfer), +// ZERO (main sequence core mass set to zero, no treatment) +// - Added new luminosity prescription for main sequence stars from Shikauchi et al. (2024) +// - Added treatment for rejuvenation of main sequence accretors when the new prescription is used -const std::string VERSION_STRING = "03.11.00"; +const std::string VERSION_STRING = "03.12.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index e09a7ad7f..c7c5f2796 100755 --- a/src/constants.h +++ b/src/constants.h @@ -290,6 +290,11 @@ constexpr double FARMER_PPISN_UPP_LIM_QUAD_REGIME = 60.0; constexpr double FARMER_PPISN_UPP_LIM_INSTABILLITY = 140.0; // Maximum CO core mass to result in PI (upper edge of PISN gap) from FARMER PPISN prescription constexpr double STARTRACK_PPISN_HE_CORE_MASS = 45.0; // Helium core mass remaining following PPISN as assumed in StarTrack (Belczynski et al. 2017 https://arxiv.org/abs/1607.03116) +constexpr double Q_CNO = 9.9073E4; // Energy released per unit mass by hydrogen fusion via the CNO cycle in Lsol Myr Msol-1 + +// Initial mass of stars above which (including the limit) we allow convective core mass calculations from Shikauchi et al. (2024) +// Note that this value should always be > 0.7 Msol +constexpr double SHIKAUCHI_LOWER_MASS_LIMIT = 15.0; // logging constants @@ -421,6 +426,7 @@ 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; + // coefficients for the calculation of initial angular frequency for Chemically Homogeneous Evolution // Mandel from Butler 2018 const DBL_VECTOR CHE_Coefficients = { 5.7914E-04, -1.9196E-06, -4.0602E-07, 1.0150E-08, -9.1792E-11, 2.9051E-13 }; @@ -3733,5 +3739,25 @@ const std::vector>> LOVERIDGE_COE }; +// Coefficients for determining Main Sequence core mass +// from Shikauchi et al. (2024), https://arxiv.org/abs/2409.00460 +// Table 2 +const std::vector SHIKAUCHI_ALPHA_COEFFICIENTS = { + {0.45, -0.0557105, -0.86589929}, // 0.1*Z_Sun + {0.45, -0.06968022, -0.73688164}, // 1/3*Z_Sun + {0.45, -0.05878711, -0.84646162} // Solar metallicity Z_Sun +}; +// Table 3 +const std::vector SHIKAUCHI_FMIX_COEFFICIENTS = { + {0.86914766, -0.60815098, 37.20654856}, // 0.1*Z_Sun + {0.86269445, -0.62623353, 35.74630996}, // 1/3*Z_Sun + {0.86605495, -0.64960375, 35.57019104} // Solar metallicity Z_Sun +}; +// Table 4 +const std::vector SHIKAUCHI_L_COEFFICIENTS = { + {3.2555795, 1.84666823, -0.79986388, -0.75728099, -0.38831172, 0.08223542, 0.49543834, 0.31314176, -0.36705796, 1.72200581}, // 0.1*Z_Sun + {3.35622529, 1.96904931, -0.88894808, -0.81112488, -0.47925922, 0.09056925, 0.53094768, 0.33971972, -0.35581284, 1.65390003}, // 1/3*Z_Sun + {3.27883249, 1.79370338, -0.71413866, -0.77019351, -0.3898752, 0.07499563, 0.5920458, 0.33846556, -0.49649838, 1.71263853} // Solar metallicity Z_Sun +}; #endif // __constants_h__ diff --git a/src/typedefs.h b/src/typedefs.h index a60d38dbc..251b7bf31 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -384,6 +384,14 @@ const COMPASUnorderedMap CHE_MODE_LABEL = { { CHE_MODE::PESSIMISTIC, "PESSIMISTIC" } }; +// main sequence core mass prescription +enum class CORE_MASS_PRESCRIPTION: int { ZERO, MANDEL, SHIKAUCHI }; +const COMPASUnorderedMap CORE_MASS_PRESCRIPTION_LABEL = { + { CORE_MASS_PRESCRIPTION::ZERO, "ZERO" }, + { CORE_MASS_PRESCRIPTION::MANDEL, "MANDEL" }, + { CORE_MASS_PRESCRIPTION::SHIKAUCHI, "SHIKAUCHI" } +}; + // logfile delimiters enum class DELIMITER: int { TAB, SPACE, COMMA }; const COMPASUnorderedMap DELIMITERLabel = { // labels diff --git a/src/yaml.h b/src/yaml.h index 274288815..ffacfd476 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -299,6 +299,7 @@ namespace yaml { " --envelope-state-prescription", " --initial-mass-function", " --LBV-mass-loss-prescription", + " --main-sequence-core-mass-prescription", " --mass-loss-prescription", " --OB-mass-loss-prescription", " --RSG-mass-loss-prescription",