Skip to content

Commit

Permalink
Merge pull request #1321 from TeamCOMPAS/1303-che-fix
Browse files Browse the repository at this point in the history 
CHE Logic Fix for Issue #1303
  • Loading branch information
@ilyamandel
ilyamandel authored Jan 14, 2025
2 parents 1676060 + 5b4f11f commit e512bd8
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 46 deletions.
104 changes: 59 additions & 45 deletions src/BaseBinaryStar.cpp
View file
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
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__

0 comments on commit e512bd8

Please sign in to comment.