From f2144a236f48b65e7b8f878616f7996b296b17e0 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Fri, 10 Jan 2025 12:26:30 -0500 Subject: [PATCH 1/6] Bug fix in CHE logic to determine whether a star is CHE based only on orbital frequency, and then set the stellar spin if required --- src/BaseBinaryStar.cpp | 4 ++-- src/changelog.h | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 3ffbfeb96..fbfd2db5e 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -310,7 +310,7 @@ void BaseBinaryStar::SetRemainingValues() { // newly-assigned rotational frequencies // star 1 - if (utils::Compare(m_Star1->Omega(), m_Star1->OmegaCHE()) >= 0) { // star 1 CH? + if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) { // star 1 CH? if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous m_Star1->SetOmega(omega); } @@ -322,7 +322,7 @@ void BaseBinaryStar::SetRemainingValues() { } // star 2 - if (utils::Compare(m_Star2->Omega(), m_Star2->OmegaCHE()) >= 0) { // star 2 CH? + if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) { // star 2 CH? if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous m_Star2->SetOmega(omega); } diff --git a/src/changelog.h b/src/changelog.h index e9c911028..c0bb5ef32 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1427,7 +1427,11 @@ // - fix for issue #1317 - SN events not always logged in BSE SN file when evolving MS merger products // - added code to ensure final BSE detailed output file TIMESTEP_COMPLETED record is always logged // (may duplicate FINAL_STATE record, but logging TIMESTEP_COMPLETED is consistent, and it's what most people look for) +// 03.10.06 VK - Jan 08, 2025 - Defect repair: +// - Fix for issue #1303 - Reduction in production of BHBH from CHE +// - Modified CHE logic to compare against orbital angular frequency rather than the stellar angular frequency, in case the stars have not been spun up yet -const std::string VERSION_STRING = "03.10.05"; + +const std::string VERSION_STRING = "03.10.06"; # endif // __changelog_h__ From 3f56742b14d291241a717dd93bd75f0f54370cde Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Fri, 10 Jan 2025 15:36:55 -0500 Subject: [PATCH 2/6] moved all CHE spin code to ProcessTides(), allowed for one or both stars to be synchronized while conserving total angular momentum --- src/BaseBinaryStar.cpp | 106 ++++++++++++++++++++++++++++++++--------- src/changelog.h | 7 +-- 2 files changed, 88 insertions(+), 25 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index fbfd2db5e..fc5a4197b 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -310,9 +310,10 @@ void BaseBinaryStar::SetRemainingValues() { // newly-assigned rotational frequencies // star 1 - if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) { // star 1 CH? - if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star1->SetOmega(omega); + if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) { // star 1 CH? + if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous + m_Star1->SetOmega(omega); } // set omega at onset of CHE } else if (m_Star1->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.7 Msol - switch if necessary @@ -322,9 +323,10 @@ void BaseBinaryStar::SetRemainingValues() { } // star 2 - if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) { // star 2 CH? - if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star2->SetOmega(omega); + if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) { // star 2 CH? + if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous + m_Star2->SetOmega(omega); } // set omega at onset of CHE } else if (m_Star2->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.0 Msol - switch if necessary @@ -1685,13 +1687,6 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_CEE); // yes - print (log) detailed output } - // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled - if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } - 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, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) @@ -2507,12 +2502,6 @@ void BaseBinaryStar::ResolveMassChanges() { m_Star2->ResetEnvelopeExpulsationByPulsations(); } - // if CHE enabled, update rotational frequency for constituent stars - assume tidally locked if tidal mode is none (otherwise, let tides do the work) - if(OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations @@ -2541,9 +2530,80 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { switch (OPTIONS->TidesPrescription()) { // which tides prescription? case TIDES_PRESCRIPTION::NONE: { // NONE - tides not enabled - // do nothing, except for CHE stars which are allowed to remain CHE even though this does not preserve angular momentum - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); + // do nothing, except for CHE stars which are allowed to remain CHE + + // if both stars are CHE, synchronize both spins and assume circular orbit, as in PERFECT tides + if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { + double omega = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), m_Star1->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU(), m_TotalAngularMomentum, OrbitalAngularVelocity()); + + if (omega >= 0.0) { // root found? (don't use utils::Compare() here) + // yes + m_Star1->SetOmega(omega); // synchronise star 1 + m_Star2->SetOmega(omega); // synchronise star 2 + + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + } + else { // no (real) root found + + // no real root found - push the binary to a common envelope + // place the constituent star closest to RLOF at RLOF and use that to + // calculate semi-major axis, then use that to calculate omega + + double ratio1 = m_Star1->StarToRocheLobeRadiusRatio(m_SemiMajorAxis, m_Star1->Mass()); // star 1 ratio of radius to Roche lobe radius + double ratio2 = m_Star2->StarToRocheLobeRadiusRatio(m_SemiMajorAxis, m_Star2->Mass()); // star 2 ratio of radius to Roche lobe radius + + double radius; + double mass1; + double mass2; + if (ratio1 >= ratio2) { // star 1 closer to RLOF than star 2 (or same)? + radius = m_Star1->Radius(); // yes - use star 1 to calculate semi-major axis at RLOF + mass1 = m_Star1->Mass(); + mass2 = m_Star2->Mass(); + } + else { // no - star 2 closer to RLOF than star 1 + radius = m_Star2->Radius(); // use star 2 to calculate semi-major axis at RLOF + mass1 = m_Star2->Mass(); + mass2 = m_Star1->Mass(); + } + + m_Eccentricity = 0.0; // assume circular + m_SemiMajorAxis = radius * RSOL_TO_AU / CalculateRocheLobeRadius_Static(mass1, mass2); // new semi-major axis - should tip into CE + } + } + // if only one star is CHE, then circularize the binary, synchronize only the CHE star, and leave the other star untouched + else if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasOneOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // one CHE star? + double che_mass; + double companion_mass; + double che_moment_of_inertia; + double Ltot; + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + che_mass = m_Star1->Mass(); + companion_mass = m_Star2->Mass(); + che_moment_of_inertia = m_Star1->CalculateMomentOfInertiaAU(); + Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity) + m_Star1->AngularMomentum(); + } + else { + che_mass = m_Star2->Mass(); + companion_mass = m_Star1->Mass(); + che_moment_of_inertia = m_Star2->CalculateMomentOfInertiaAU(); + Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity) + m_Star2->AngularMomentum(); + } + // Use the synchronization formula assuming a circular orbit and ignoring the companion angular momentum + double omega = OmegaAfterSynchronisation(che_mass, companion_mass, che_moment_of_inertia, 0.0, Ltot, OrbitalAngularVelocity()); + if (omega >= 0.0) { // root found? (don't use utils::Compare() here) + // yes + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} + else if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} + + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + } + } } break; @@ -2569,6 +2629,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { m_SemiMajorAxis = m_SemiMajorAxis + ((DSemiMajorAxis1Dt + DSemiMajorAxis2Dt) * p_Dt * MYR_TO_YEAR); // evolve separation m_Eccentricity = m_Eccentricity + ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR); // evolve eccentricity m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } break; @@ -2587,6 +2648,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis m_Eccentricity = 0.0; // circularise m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } else { // no (real) root found diff --git a/src/changelog.h b/src/changelog.h index c0bb5ef32..fdb02f949 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1427,11 +1427,12 @@ // - fix for issue #1317 - SN events not always logged in BSE SN file when evolving MS merger products // - added code to ensure final BSE detailed output file TIMESTEP_COMPLETED record is always logged // (may duplicate FINAL_STATE record, but logging TIMESTEP_COMPLETED is consistent, and it's what most people look for) -// 03.10.06 VK - Jan 08, 2025 - Defect repair: +// 03.11.00 VK - Jan 08, 2025 - Enhancement, Defect repair: // - Fix for issue #1303 - Reduction in production of BHBH from CHE -// - Modified CHE logic to compare against orbital angular frequency rather than the stellar angular frequency, in case the stars have not been spun up yet +// - Modified CHE initialization logic to compare against orbital angular frequency rather than the stellar angular frequency, in case the stars have not been spun up yet +// - Moved all CHE rotation related code to ProcessTides(), ensuring that any spin up conserves total angular momentum -const std::string VERSION_STRING = "03.10.06"; +const std::string VERSION_STRING = "03.11.00"; # endif // __changelog_h__ From 50d6c9954751f0fd15b6d50cffb9ffb1e3d189d6 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Fri, 10 Jan 2025 15:56:17 -0500 Subject: [PATCH 3/6] Optimized CHE treatment to handle one or more CHE stars with same code --- src/BaseBinaryStar.cpp | 79 ++++++++++-------------------------------- 1 file changed, 19 insertions(+), 60 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index fc5a4197b..67897baa3 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2532,77 +2532,36 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { // do nothing, except for CHE stars which are allowed to remain CHE - // if both stars are CHE, synchronize both spins and assume circular orbit, as in PERFECT tides - if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasTwoOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { - double omega = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), m_Star1->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU(), m_TotalAngularMomentum, OrbitalAngularVelocity()); + // if at least one star is CHE, then circularize the binary and synchronize only the CHE stars conserving total angular momentum + if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasOneOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // one CHE star? + double che_I1 = 0.0; + double che_I2 = 0.0; + double che_Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); - if (omega >= 0.0) { // root found? (don't use utils::Compare() here) - // yes - m_Star1->SetOmega(omega); // synchronise star 1 - m_Star2->SetOmega(omega); // synchronise star 2 - - m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis - m_Eccentricity = 0.0; // circularise - m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum - m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); - } - else { // no (real) root found - - // no real root found - push the binary to a common envelope - // place the constituent star closest to RLOF at RLOF and use that to - // calculate semi-major axis, then use that to calculate omega - - double ratio1 = m_Star1->StarToRocheLobeRadiusRatio(m_SemiMajorAxis, m_Star1->Mass()); // star 1 ratio of radius to Roche lobe radius - double ratio2 = m_Star2->StarToRocheLobeRadiusRatio(m_SemiMajorAxis, m_Star2->Mass()); // star 2 ratio of radius to Roche lobe radius - - double radius; - double mass1; - double mass2; - if (ratio1 >= ratio2) { // star 1 closer to RLOF than star 2 (or same)? - radius = m_Star1->Radius(); // yes - use star 1 to calculate semi-major axis at RLOF - mass1 = m_Star1->Mass(); - mass2 = m_Star2->Mass(); - } - else { // no - star 2 closer to RLOF than star 1 - radius = m_Star2->Radius(); // use star 2 to calculate semi-major axis at RLOF - mass1 = m_Star2->Mass(); - mass2 = m_Star1->Mass(); - } - - m_Eccentricity = 0.0; // assume circular - m_SemiMajorAxis = radius * RSOL_TO_AU / CalculateRocheLobeRadius_Static(mass1, mass2); // new semi-major axis - should tip into CE - } - } - // if only one star is CHE, then circularize the binary, synchronize only the CHE star, and leave the other star untouched - else if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasOneOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // one CHE star? - double che_mass; - double companion_mass; - double che_moment_of_inertia; - double Ltot; if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { - che_mass = m_Star1->Mass(); - companion_mass = m_Star2->Mass(); - che_moment_of_inertia = m_Star1->CalculateMomentOfInertiaAU(); - Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity) + m_Star1->AngularMomentum(); + che_I1 = m_Star1->CalculateMomentOfInertiaAU(); + che_Ltot += m_Star1->AngularMomentum(); } - else { - che_mass = m_Star2->Mass(); - companion_mass = m_Star1->Mass(); - che_moment_of_inertia = m_Star2->CalculateMomentOfInertiaAU(); - Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity) + m_Star2->AngularMomentum(); + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + che_I2 = m_Star2->CalculateMomentOfInertiaAU(); + che_Ltot += m_Star2->AngularMomentum(); } - // Use the synchronization formula assuming a circular orbit and ignoring the companion angular momentum - double omega = OmegaAfterSynchronisation(che_mass, companion_mass, che_moment_of_inertia, 0.0, Ltot, OrbitalAngularVelocity()); - if (omega >= 0.0) { // root found? (don't use utils::Compare() here) - // yes + + double omega = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), che_I1, che_I2, che_Ltot, OrbitalAngularVelocity()); + if (omega >= 0.0) { // root found? (don't use utils::Compare() here) + // yes if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} - else if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis m_Eccentricity = 0.0; // circularise m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } + else { // no (real) root found, synchronize CHE stars ignoring angular momentum conservation + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} + } } } break; From 3e1066b7ff730c9472ad188187658cc01b5571c9 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Tue, 14 Jan 2025 11:58:29 -0500 Subject: [PATCH 4/6] Updated changelog to highlight CHE spin assumptions --- src/changelog.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/changelog.h b/src/changelog.h index a8611a5f4..255cc8f4d 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1430,9 +1430,11 @@ // 03.10.06 VK - Jan 13, 2025 - Enhancement: // - Modified the KAPIL2024 tides to ignore quadratic 'e' terms (for spin and separation evolution) if they spin up an already synchronized star. // 03.11.00 VK - Jan 14, 2025 - Enhancement, Defect repair: -// - Fix for issue #1303 - Reduction in production of BHBH from CHE -// - Modified CHE initialization logic to compare against orbital angular frequency rather than the stellar angular frequency, in case the stars have not been spun up yet -// - Moved all CHE rotation related code to ProcessTides(), ensuring that any spin up conserves total angular momentum +// - Fix for issue #1303 - Reduction in production of BHBH from CHE, other CHE-related improvements. +// - Stars that have sufficiently rapid angular frequencies at ZAMS are now initialized as CHE stars, regardless of the tidal prescription. +// - 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. const std::string VERSION_STRING = "03.11.00"; From 4c77d2a77f8d8de70e2c1f8e33646e43604ff418 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Tue, 14 Jan 2025 11:59:47 -0500 Subject: [PATCH 5/6] if stellar rotation is user specified, do not modify spin at CHE initialization --- src/BaseBinaryStar.cpp | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 0acfe84d9..d99b68d92 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -301,19 +301,14 @@ void BaseBinaryStar::SetRemainingValues() { // check for CHE // - // because we've changed the rotational frequency of the constituent stars we - // have to reset the stellar type - at this stage, based on their rotational - // frequency at birth, they may have already been assigned one of MS_LTE_07, - // MS_GT_07, or CHEMICALLY_HOMOGENEOUS - // - // here we need to change from MS_* -> CH, or from CH->MS* based on the - // newly-assigned rotational frequencies + // here we need to change from MS_* -> CH, or from CH->MS* based on the binary orbital frequency, assuming that the stars are tidally locked + // set the spin to the orbital frequency, unless the user has specified a spin frequency // star 1 if (utils::Compare(omega, m_Star1->OmegaCHE()) >= 0) { // star 1 CH? if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star1->SetOmega(omega); } // set omega at onset of CHE + if (OPTIONS->OptionDefaulted("rotational-frequency-1")) m_Star1->SetOmega(omega); } // set spin to orbital frequency unless user specified } else if (m_Star1->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.7 Msol - switch if necessary @@ -326,7 +321,7 @@ void BaseBinaryStar::SetRemainingValues() { if (utils::Compare(omega, m_Star2->OmegaCHE()) >= 0) { // star 2 CH? if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true); // yes, switch if not already Chemically Homogeneous - m_Star2->SetOmega(omega); } // set omega at onset of CHE + if (OPTIONS->OptionDefaulted("rotational-frequency-2")) m_Star2->SetOmega(omega); } // set spin to orbital frequency unless user specified } else if (m_Star2->MZAMS() <= 0.7) { // no - MS - initial mass determines actual type (don't use utils::Compare() here) if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true); // MS <= 0.0 Msol - switch if necessary @@ -2558,13 +2553,13 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { che_Ltot += m_Star2->AngularMomentum(); } - double omega = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), che_I1, che_I2, che_Ltot, OrbitalAngularVelocity()); - if (omega >= 0.0) { // root found? (don't use utils::Compare() here) + double omega_sync = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), che_I1, che_I2, che_Ltot, omega); + if (omega_sync >= 0.0) { // root found? (don't use utils::Compare() here) // yes - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega_sync);} + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega_sync);} - m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega_sync / omega_sync); // re-calculate semi-major axis m_Eccentricity = 0.0; // circularise m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); From 5b4f11fc9a3c6327c0acf3f2636f17af37f509ec Mon Sep 17 00:00:00 2001 From: jeffriley Date: Wed, 15 Jan 2025 09:39:21 +1100 Subject: [PATCH 6/6] Cleanup --- src/BaseBinaryStar.cpp | 54 ++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 28 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index d99b68d92..a451b0658 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2396,7 +2396,7 @@ double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis, double Jorb = CalculateOrbitalAngularMomentum(p_Star1Mass, p_Star2Mass, p_SemiMajorAxis, p_Eccentricity); - return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb; + return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb; } @@ -2418,8 +2418,8 @@ void BaseBinaryStar::CalculateEnergyAndAngularMomentum() { m_OrbitalAngularMomentumPrev = m_OrbitalAngularMomentum; m_TotalAngularMomentumPrev = m_TotalAngularMomentum; - double totalMass = m_Star1->Mass() + m_Star2->Mass(); - double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; + double totalMass = m_Star1->Mass() + m_Star2->Mass(); + double reducedMass = (m_Star1->Mass() * m_Star2->Mass()) / totalMass; m_OrbitalEnergy = CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis); m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); @@ -2457,9 +2457,9 @@ void BaseBinaryStar::ResolveMassChanges() { if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius - double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? - massChange * sqrt(G_AU_Msol_yr * m_Star1->Mass() * m_Star1->Radius() * RSOL_TO_AU) : - (2.0 / 3.0) * massChange * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Omega(); + double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) + ? massChange * sqrt(G_AU_Msol_yr * m_Star1->Mass() * m_Star1->Radius() * RSOL_TO_AU) + : (2.0 / 3.0) * massChange * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Radius() * RSOL_TO_AU * m_Star1->Omega(); // update mass of star according to mass loss and mass transfer, then update age accordingly (void)m_Star1->UpdateAttributes(massChange, 0.0); // update mass for star m_Star1->UpdateInitialMass(); // update effective initial mass of star (MS, HG & HeMS) @@ -2480,9 +2480,9 @@ void BaseBinaryStar::ResolveMassChanges() { double massChange = m_Star2->MassLossDiff() + m_Star2->MassTransferDiff(); // mass change due to winds and mass transfer if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius - double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? - massChange * sqrt(G_AU_Msol_yr * m_Star2->Mass() * m_Star2->Radius() * RSOL_TO_AU) : - (2.0 / 3.0) * massChange * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Omega(); + double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) + ? massChange * sqrt(G_AU_Msol_yr * m_Star2->Mass() * m_Star2->Radius() * RSOL_TO_AU) + : (2.0 / 3.0) * massChange * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Radius() * RSOL_TO_AU * m_Star2->Omega(); // update mass of star according to mass loss and mass transfer, then update age accordingly (void)m_Star2->UpdateAttributes(massChange, 0.0); // update mass for star m_Star2->UpdateInitialMass(); // update effective initial mass of star (MS, HG & HeMS) @@ -2508,8 +2508,6 @@ void BaseBinaryStar::ResolveMassChanges() { m_Star2->ResetEnvelopeExpulsationByPulsations(); } - - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations if ((m_Star1->StellarType() != stellarType1) || (m_Star2->StellarType() != stellarType2)) { // stellar type change? @@ -2530,7 +2528,6 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { if (!m_Unbound) { // binary bound? // yes - process tides if enabled - double omega = OrbitalAngularVelocity(); switch (OPTIONS->TidesPrescription()) { // which tides prescription? @@ -2540,31 +2537,32 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { // if at least one star is CHE, then circularize the binary and synchronize only the CHE stars conserving total angular momentum if (OPTIONS->CHEMode() != CHE_MODE::NONE && HasOneOf({STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS})) { // one CHE star? - double che_I1 = 0.0; - double che_I2 = 0.0; + double che_I1 = 0.0; + double che_I2 = 0.0; double che_Ltot = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { - che_I1 = m_Star1->CalculateMomentOfInertiaAU(); + che_I1 = m_Star1->CalculateMomentOfInertiaAU(); che_Ltot += m_Star1->AngularMomentum(); - } + } + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { - che_I2 = m_Star2->CalculateMomentOfInertiaAU(); + che_I2 = m_Star2->CalculateMomentOfInertiaAU(); che_Ltot += m_Star2->AngularMomentum(); - } + } double omega_sync = OmegaAfterSynchronisation(m_Star1->Mass(), m_Star2->Mass(), che_I1, che_I2, che_Ltot, omega); - if (omega_sync >= 0.0) { // root found? (don't use utils::Compare() here) - // yes + if (omega_sync >= 0.0) { // root found? (don't use utils::Compare() here) + // yes if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega_sync);} if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega_sync);} - m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega_sync / omega_sync); // re-calculate semi-major axis - m_Eccentricity = 0.0; // circularise - m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega_sync / omega_sync); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } - else { // no (real) root found, synchronize CHE stars ignoring angular momentum conservation + else { // no (real) root found, synchronize CHE stars ignoring angular momentum conservation if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star1->SetOmega(omega);} if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS){m_Star2->SetOmega(omega);} } @@ -2610,9 +2608,9 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { m_Star1->SetOmega(omega); // synchronise star 1 m_Star2->SetOmega(omega); // synchronise star 2 - m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis - m_Eccentricity = 0.0; // circularise - m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum + m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / omega / omega); // re-calculate semi-major axis + m_Eccentricity = 0.0; // circularise + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); } else { // no (real) root found @@ -2656,6 +2654,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { } } + /* * Root solver to determine rotational frequency after synchronisation for tides * @@ -2759,7 +2758,6 @@ double BaseBinaryStar::OmegaAfterSynchronisation(const double p_M1, const double } - /* * Calculate and emit gravitational radiation. *