Skip to content
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

CHE Logic Fix for Issue #1303 #1321

Merged
merged 7 commits into from
Jan 14, 2025
Merged

CHE Logic Fix for Issue #1303 #1321

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f2144a2
Bug fix in CHE logic to determine whether a star is CHE based only on…
veome22 Jan 10, 2025
3f56742
moved all CHE spin code to ProcessTides(), allowed for one or both st…
veome22 Jan 10, 2025
50d6c99
Optimized CHE treatment to handle one or more CHE stars with same code
veome22 Jan 10, 2025
4fb9d55
Merge branch 'dev' into 1303-che-fix
veome22 Jan 14, 2025
3e1066b
Updated changelog to highlight CHE spin assumptions
veome22 Jan 14, 2025
4c77d2a
if stellar rotation is user specified, do not modify spin at CHE init…
veome22 Jan 14, 2025
5b4f11f
Cleanup
jeffriley Jan 14, 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
104 changes: 59 additions & 45 deletions src/BaseBinaryStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;
}


Expand All @@ -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);
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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?
Expand All @@ -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;

Expand All @@ -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;

Expand All @@ -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

Expand Down Expand Up @@ -2640,6 +2654,7 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) {
}
}


/*
* Root solver to determine rotational frequency after synchronisation for tides
*
Expand Down Expand Up @@ -2743,7 +2758,6 @@ double BaseBinaryStar::OmegaAfterSynchronisation(const double p_M1, const double
}



/*
* Calculate and emit gravitational radiation.
*
Expand Down
8 changes: 7 additions & 1 deletion src/changelog.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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__
Loading