Skip to content

Fixes to CHE and MS mergers with BRCEK core mass prescription #1401

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
brcekadam merged 6 commits into from
Jun 23, 2025
Merged

Fixes to CHE and MS mergers with BRCEK core mass prescription #1401

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a07ae2a
Fixes and upgrades to CHE and MS mergers
brcekadam Jun 20, 2025
880ca51
Update changelog.h
brcekadam Jun 23, 2025
f8b1be0
Typo fix
brcekadam Jun 23, 2025
c23bc48
Fix to hydrogen mass calculations
brcekadam Jun 23, 2025
02468cd
Typo fix
brcekadam Jun 23, 2025
5ab3409
Fix to hydrogen mass calculation
brcekadam Jun 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 6 additions & 7 deletions src/BaseBinaryStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1778,15 +1778,14 @@ void BaseBinaryStar::ResolveMainSequenceMerger() {
(utils::Compare(m_Star1->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0) &&
(utils::Compare(m_Star2->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)) {

double coreMass1 = m_Star1->MainSequenceCoreMass();
double coreMass2 = m_Star2->MainSequenceCoreMass();
// total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope)
double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + m_Star1->HydrogenAbundanceSurface() * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0);
double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + m_Star2->HydrogenAbundanceSurface() * std::max(m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass(), 0.0);

double coreHeliumMass1 = m_Star1->HeliumAbundanceCore() * coreMass1;
double coreHeliumMass2 = m_Star2->HeliumAbundanceCore() * coreMass2;

finalHydrogenMass = finalMass * initialHydrogenFraction - coreHeliumMass1 - coreHeliumMass2;
// final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope of the star with higher surface hydrogen abundance
finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * std::max(m_Star1->HydrogenAbundanceSurface(), m_Star2->HydrogenAbundanceSurface());
}
else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction;
else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction;

m_SemiMajorAxis = std::numeric_limits<float>::infinity(); // set separation to infinity to avoid subsequent fake interactions with a massless companion (RLOF, CE, etc.)

Expand Down
1 change: 1 addition & 0 deletions src/BaseStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
// initialise remaining member variables

// Zero age main sequence parameters
m_InitialMainSequenceCoreMass = DEFAULT_INITIAL_DOUBLE_VALUE; // initialised in MS_gt_07 class if BRCEK core mass prescription is used
m_RZAMS = CalculateRadiusAtZAMS(m_MZAMS);
m_LZAMS = CalculateLuminosityAtZAMS(m_MZAMS);
m_TZAMS = CalculateTemperatureOnPhase_Static(m_LZAMS, m_RZAMS);
Expand Down
2 changes: 2 additions & 0 deletions src/BaseStar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,7 @@ class BaseStar {
double HydrogenAbundanceSurface() const { return m_HydrogenAbundanceSurface; }
double InitialHeliumAbundance() const { return m_InitialHeliumAbundance; }
double InitialHydrogenAbundance() const { return m_InitialHydrogenAbundance; }
double InitialMainSequenceCoreMass() const { return m_InitialMainSequenceCoreMass; }
bool IsAIC() const { return (m_SupernovaDetails.events.current & SN_EVENT::AIC) == SN_EVENT::AIC; }
bool IsCCSN() const { return (m_SupernovaDetails.events.current & SN_EVENT::CCSN) == SN_EVENT::CCSN; }
bool IsHeSD() const { return (m_SupernovaDetails.events.current & SN_EVENT::HeSD) == SN_EVENT::HeSD; }
Expand Down Expand Up @@ -404,6 +405,7 @@ class BaseStar {
// Zero Age Main Sequence
double m_InitialHeliumAbundance; // Initial helium abundance (Y)
double m_InitialHydrogenAbundance; // Initial hydrogen abundance (X)
double m_InitialMainSequenceCoreMass; // Initial main sequence core mass (used in BRCEK core mass prescription)
double m_LZAMS; // ZAMS Luminosity
double m_MZAMS; // ZAMS Mass
double m_OmegaZAMS; // ZAMS Angular Frequency
Expand Down
4 changes: 4 additions & 0 deletions src/CH.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -423,6 +423,10 @@ STELLAR_TYPE CH::EvolveToNextPhase() {
if (m_Age < m_Timescales[static_cast<int>(TIMESCALE::tMS)]) { // evolving off because of age?
stellarType = STELLAR_TYPE::MS_GT_07; // no - must have spun down - evolve as MS star now
m_CHE = false; // evolved CH->MS

// if BRCEK core mass calculations enabled, initialise the core mass based on current mass and central helium fraction
if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0))
m_MainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_Mass, m_HeliumAbundanceCore);
}
else { // yes
stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_MS; // evolve as HeMS star now
Expand Down
2 changes: 2 additions & 0 deletions src/CH.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,8 @@ class CH: virtual public BaseStar, public MS_gt_07 {
bool ShouldEvolveOnPhase() const { return m_Age < m_Timescales[static_cast<int>(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 UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { }; // Do not use core mass calculations during CHE phase

};

Expand Down
2 changes: 1 addition & 1 deletion src/MS_gt_07.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ class MS_gt_07: virtual public BaseStar, public MainSequence {
utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0 && // ZAMS mass >= BRCEK_LOWER_MASS_LIMIT?
m_Time <= 0.0) { // star not yet aged past creation?
// yes - initialise
m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS);
m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS, m_InitialHeliumAbundance);
m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass;
m_Luminosity = MainSequence::CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0);
m_Radius = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Tau, m_RZAMS0);
Expand Down
41 changes: 28 additions & 13 deletions src/MainSequence.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -289,8 +289,9 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function

// If BRCEK core prescription is used, return luminosity from Shikauchi et al. (2024) during core hydrogen burning (valid for MZAMS >= 15 Msol) or
// luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT)
if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) {
// luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT); do not use Shikauchi luminosity
// prescription during CHE
if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0) && !m_CHE) {
if (utils::Compare(p_Time, 0.99 * timescales(tMS)) > 0) // star in MS hook?
return CalculateLuminosityTransitionToHG(p_Mass, p_Time, p_LZAMS);
else {
Expand Down Expand Up @@ -826,18 +827,30 @@ DBL_DBL MainSequence::CalculateMainSequenceCoreMassBrcek(const double p_Dt, cons


/*
* 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
* Calculate the initial convective core mass of a main sequence star at ZAMS using Equation (A3) from Shikauchi+ (2024),
* or after full mixing (due to merger or CHE) for an arbitrary central helium fraction using the approach
* described in Brcek+ (2025)
*
* double CalculateInitialMainSequenceCoreMass(const double p_MZAMS)
* double CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore)
*
* @param [IN] p_MZAMS Mass at ZAMS or after merger in Msol
* @param [IN] p_Mass Mass at ZAMS, after merger or after spin down of CH star in Msol
* @param [IN] p_HeliumAbundanceCore Central helium fraction
* @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;
double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) const {

double fmix = 0.0;
// At ZAMS, use the approach from Shikauchi+ (2024)
if (utils::Compare(p_HeliumAbundanceCore, m_InitialHeliumAbundance) == 0) {
DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS);
fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_Mass / fmixCoefficients[2]);
}
// After full mixing not at ZAMS, use the approach from Brcek+ (2025)
else {
double h = PPOW(10.0, p_HeliumAbundanceCore * (p_HeliumAbundanceCore + 2.0) / 4.0);
fmix = BRCEK_FMIX_COEFFICIENTS[0] + BRCEK_FMIX_COEFFICIENTS[1] * std::exp(-p_Mass * h / BRCEK_FMIX_COEFFICIENTS[2]) * PPOW(1.0 - BRCEK_FMIX_COEFFICIENTS[4] / (p_Mass * h), BRCEK_FMIX_COEFFICIENTS[3]);
}
return fmix * p_Mass;
}


Expand Down Expand Up @@ -1142,12 +1155,14 @@ void MainSequence::UpdateAfterMerger(double p_Mass, double p_HydrogenMass) {
m_Age = m_Tau * timescales(tMS);

m_HeliumAbundanceCore = 1.0 - m_Metallicity - p_HydrogenMass / p_Mass;

m_HydrogenAbundanceCore = 1.0 - m_Metallicity - m_HeliumAbundanceCore;

m_HeliumAbundanceSurface = m_HeliumAbundanceCore; // abundances are the same throughout the star, assuming uniform mixing after merger
m_HydrogenAbundanceSurface = m_HydrogenAbundanceCore;

if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) {
m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass); // update initial mixing core mass
m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; // update core mass
m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass, m_HeliumAbundanceCore); // update initial mixing core mass
m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; // update core mass
}

UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true);
Expand Down
4 changes: 1 addition & 3 deletions src/MainSequence.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,7 @@ class MainSequence: virtual public BaseStar {
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;
Expand Down Expand Up @@ -78,7 +76,7 @@ class MainSequence: virtual public BaseStar {
double CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore) const;
double CalculateLuminosityTransitionToHG(const double p_Mass, const double p_Age, double const p_LZAMS) const;
DBL_DBL CalculateMainSequenceCoreMassBrcek(const double p_Dt, const double p_MassLossRate);
double CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const;
double CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) 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)
Expand Down
1 change: 1 addition & 0 deletions src/Star.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ class Star {
double HydrogenAbundanceSurface() const { return m_Star->HydrogenAbundanceSurface(); }
double InitialHeliumAbundance() const { return m_Star->InitialHeliumAbundance(); }
double InitialHydrogenAbundance() const { return m_Star->InitialHydrogenAbundance(); }
double InitialMainSequenceCoreMass() const { return m_Star->InitialMainSequenceCoreMass(); }
bool IsAIC() const { return m_Star->IsAIC(); }
bool IsCCSN() const { return m_Star->IsCCSN(); }
bool IsDegenerate() const { return m_Star->IsDegenerate(); }
Expand Down
5 changes: 4 additions & 1 deletion src/changelog.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1583,8 +1583,11 @@
// - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02)
// - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static()
// - Removed McSN from GBParams, instead computed on the fly when needed
// 03.20.04 AB - Jun 23, 2025 - Defect repair, enhancement:
// - Fixes to MS mergers and CHE when BRCEK core mass prescription is used -- MS core mass is now correctly initialised after full mixing in MS
// mergers and CH stars that spun down

const std::string VERSION_STRING = "03.20.03";
const std::string VERSION_STRING = "03.20.04";


# endif // __changelog_h__
4 changes: 4 additions & 0 deletions src/constants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3811,5 +3811,9 @@ const std::vector<DBL_VECTOR> SHIKAUCHI_L_COEFFICIENTS = {
{3.45464814, 0.94880846, -1.11409154, -0.86672079, 2.38986855, 0.04448855, -2.74913945, 0.60905625, -0.27648361, 0.03514139, 1.37569819, -0.19184532, 0.12816567, -0.14392935, 1.76390159}, // 1/3*Z_Sun
{3.80166901, 0.37407948, -1.29904749, -1.34541622, 3.70934166, 0.28320469, -3.92327169, 0.92444477, -0.40146717, -0.00821364, 1.80297947, -0.15776603, 0.09205681, -0.21913557, 1.78496679} // Solar metallicity Z_Sun
};
// Coefficients used to determine the initial convective core mass of MS star after full mixing (due to a merger or CHE)
// from Brcek et al. (2025)
const DBL_VECTOR BRCEK_FMIX_COEFFICIENTS = {0.897487516925644, -0.570499075062176, 56.9817940103868, 0.54483109433479, 1.72842835324164};


#endif // __constants_h__