diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index e405d8f35..a451b0658 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -301,18 +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(m_Star1->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 + 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 @@ -322,9 +318,10 @@ void BaseBinaryStar::SetRemainingValues() { } // star 2 - if (utils::Compare(m_Star2->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 + 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 @@ -1696,13 +1693,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) @@ -2406,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; } @@ -2428,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); @@ -2467,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) @@ -2490,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) @@ -2518,14 +2508,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 if ((m_Star1->StellarType() != stellarType1) || (m_Star2->StellarType() != stellarType2)) { // stellar type change? @@ -2546,15 +2528,45 @@ 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? 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 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 (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + che_I1 = m_Star1->CalculateMomentOfInertiaAU(); + che_Ltot += m_Star1->AngularMomentum(); + } + + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) { + 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 (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_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; @@ -2580,6 +2592,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; @@ -2595,9 +2608,10 @@ 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 @@ -2640,6 +2654,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { } } + /* * Root solver to determine rotational frequency after synchronisation for tides * @@ -2743,7 +2758,6 @@ double BaseBinaryStar::OmegaAfterSynchronisation(const double p_M1, const double } - /* * Calculate and emit gravitational radiation. * diff --git a/src/changelog.h b/src/changelog.h index 409f94b7e..255cc8f4d 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1429,7 +1429,13 @@ // (may duplicate FINAL_STATE record, but logging TIMESTEP_COMPLETED is consistent, and it's what most people look for) // 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, 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.10.06"; +const std::string VERSION_STRING = "03.11.00"; # endif // __changelog_h__