From 1a542bae25c90e6bd7a35383f4937fdc25862723 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Tue, 7 Jan 2025 15:20:04 -0500 Subject: [PATCH 1/3] Modified tides to ignore terms quadratic in eccentricity, when they push a star to spin beyond synchronously --- src/BaseBinaryStar.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 93f4a033e..587fd4802 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1052,7 +1052,8 @@ double BaseBinaryStar::CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImK double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_8 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - + + // No need to ignore quadratic e order terms during (super) synchronous rotation, since this formula is already linear in eccentricity return -(3.0 / 4.0) * (m_Eccentricity / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * ((3.0 * ImK10 / 2.0) - (ImK12 / 4.0) - ImK22 + (49.0 * ImK32 / 4.0)); } @@ -1080,8 +1081,11 @@ double BaseBinaryStar::CalculateDOmegaTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, con double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_6 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; + double e2_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); - return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + ((m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)))); + // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_term), 0.0) > 0)){e2_term = 0.0;} + return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + e2_term); } @@ -1108,8 +1112,13 @@ double BaseBinaryStar::CalculateDSemiMajorAxisTidalDt(const DBL_DBL_DBL_DBL p_Im double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_7 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; + double e2_term = (m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)); + + // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms + double e2_term_spin = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_term_spin), 0.0) > 0)){e2_term = 0.0;} - return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + ((m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)))); + return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + e2_term); } From fe9cf1e112b7538a5810ac07f7923b1551dcadce Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Tue, 7 Jan 2025 15:27:01 -0500 Subject: [PATCH 2/3] updated changelog --- src/changelog.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/changelog.h b/src/changelog.h index fc9d05f50..effe48016 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1423,7 +1423,9 @@ // - fix for issue #1310 - run terminates prematurely if error in grid file // 03.10.04 RTW - Nov 27, 2024 - Defect repair: // - fix for issue #1247 - SN Euler angles had incomplete logic, leading to a div by zero in some cases +// 03.10.05 VK - Jan 7, 2025 - Enhancement: +// - Modified the KAPIL2024 tides to ignore quadratic 'e' terms for spin and separation evolution if they incorrectly spin up an already synchronized star. -const std::string VERSION_STRING = "03.10.04"; +const std::string VERSION_STRING = "03.10.05"; # endif // __changelog_h__ From c5d1d7ff04fb40a8e06be2cf22129193bb388b4c Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Mon, 13 Jan 2025 20:50:07 -0500 Subject: [PATCH 3/3] clarified when the higher order terms should be ignored in the comments --- src/BaseBinaryStar.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ebc5badeb..e405d8f35 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1081,11 +1081,11 @@ double BaseBinaryStar::CalculateDOmegaTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, con double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_6 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - double e2_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + double e2_spin_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms - if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_term), 0.0) > 0)){e2_term = 0.0;} - return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + e2_term); + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_spin_term), 0.0) > 0)){e2_spin_term = 0.0;} + return (3.0 / 2.0) * (1.0 / MoIstar) * (G_AU_Msol_yr * massCompanion * massCompanion / R1_AU) * R1_over_a_6 * (ImK22 + e2_spin_term); } @@ -1112,13 +1112,15 @@ double BaseBinaryStar::CalculateDSemiMajorAxisTidalDt(const DBL_DBL_DBL_DBL p_Im double R1_AU = radiusStar * RSOL_TO_AU; double R1_over_a = R1_AU / m_SemiMajorAxis; double R1_over_a_7 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - double e2_term = (m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)); + double e2_sma_term = (m_Eccentricity * m_Eccentricity) * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)); - // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms - double e2_term_spin = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); - if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_term_spin), 0.0) > 0)){e2_term = 0.0;} + // if the star is rotating (super) synchronously AND quadratic 'e' terms cause the star to spin up further, ignore the higher order terms. + // Note: here we use the SPIN e^2 terms (not the semi-major axis terms) to determine when to ignore the higher order terms in semi-major axis evolution. + // this is to ensure that the higher order terms are always consistently applied/ignored across the tidal evolution equations. + double e2_spin_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + if ((utils::Compare(p_Star->Omega(), OrbitalAngularVelocity()) > 0) && (utils::Compare((ImK22 + e2_spin_term), 0.0) > 0)){e2_sma_term = 0.0;} - return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + e2_term); + return -(3.0 / OrbitalAngularVelocity()) * (1.0 + (massCompanion / massStar)) * (G_AU_Msol_yr * massCompanion / R1_AU / R1_AU) * R1_over_a_7 * (ImK22 + e2_sma_term); }