From fd038ced63f85073b38fc4da051b24529a9ce225 Mon Sep 17 00:00:00 2001 From: Veome Kapil Date: Mon, 8 Apr 2024 14:41:11 +1000 Subject: [PATCH] Generalized tidal love number and resulting secular equations to other (l,m) modes --- src/BaseBinaryStar.cpp | 22 ++++++++++++++-------- src/BaseBinaryStar.h | 27 ++++++++++++++++++--------- src/BaseStar.h | 2 +- src/GiantBranch.cpp | 20 +++++++++++--------- src/GiantBranch.h | 2 +- src/Star.h | 2 +- 6 files changed, 46 insertions(+), 29 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index c81f7511f..623a8b207 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2335,7 +2335,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations if (OPTIONS->EnableRealisticTides() && !m_Unbound) { - // Change binary semi-major axis and spin of each star based on realistic tidal torque + // Change binary semi-major axis, eccentricity, and spin of each star based on realistic tidal torque // Adjust the binary orbital frequency to match semi-major axis. // if m_Omega == 0.0 (should only happen on the first timestep), calculate m_Omega here if (utils::Compare(m_Omega, 0.0) == 0) { @@ -2343,7 +2343,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { } - double DSemiMajorAxis1Dt = CalculateDSemiMajorAxisTidalDt(m_Star1->CalculateImK22Tidal(m_Omega), + double DSemiMajorAxis1Dt = CalculateDSemiMajorAxisTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star1->Mass(), m_Star1->Radius(), m_Star2->Mass(), @@ -2351,7 +2352,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(m_Star2->CalculateImK22Tidal(m_Omega), + double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star2->Mass(), m_Star2->Radius(), m_Star1->Mass(), @@ -2359,7 +2361,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DEccentricity1Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImK22Tidal(m_Omega), + double DEccentricity1Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star1->Mass(), m_Star1->Radius(), m_Star2->Mass(), @@ -2367,7 +2370,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DEccentricity2Dt = CalculateDEccentricityTidalDt(m_Star2->CalculateImK22Tidal(m_Omega), + double DEccentricity2Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star2->Mass(), m_Star2->Radius(), m_Star1->Mass(), @@ -2375,7 +2379,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DOmega1Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImK22Tidal(m_Omega), + double DOmega1Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star1->Mass(), m_Star1->Radius(), m_Star1->CalculateMomentOfInertiaAU(), @@ -2384,7 +2389,8 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis, m_Eccentricity); - double DOmega2Dt = CalculateDOmegaTidalDt(m_Star2->CalculateImK22Tidal(m_Omega), + double DOmega2Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImKlmTidal(m_Omega, 1, 0), m_Star1->CalculateImKlmTidal(m_Omega, 1, 2), + m_Star1->CalculateImKlmTidal(m_Omega, 2, 2), m_Star1->CalculateImKlmTidal(m_Omega, 3, 2), m_Star2->Mass(), m_Star2->Radius(), m_Star2->CalculateMomentOfInertiaAU(), @@ -2398,7 +2404,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { m_SemiMajorAxis = m_SemiMajorAxis + ((DSemiMajorAxis1Dt + DSemiMajorAxis2Dt) * p_Dt * MYR_TO_YEAR); // change separation - m_Eccentricity = m_Eccentricity + ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR); // change eccentricity (This increases eccentricity for positive k22) + m_Eccentricity = m_Eccentricity + ((DEccentricity1Dt + DEccentricity2Dt) * p_Dt * MYR_TO_YEAR); // change eccentricity m_Omega = std::sqrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / m_SemiMajorAxis / m_SemiMajorAxis / m_SemiMajorAxis); // re-calculate orbital frequency diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 0bdd2eaed..7ea37d3d4 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -774,9 +774,12 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.6) * * - * double CalculateDSemiMajorAxisTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDSemiMajorAxisTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * + * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_M2 Mass of companion star (Msol) @@ -785,12 +788,12 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in semi-major axis for binary (AU/yr) */ - double CalculateDSemiMajorAxisTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDSemiMajorAxisTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_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; - return - (3.0 / p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2/ R1_AU / R1_AU) * R1_over_a_7 * p_ImK22; + return - (3.0 / p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2/ R1_AU / R1_AU) * R1_over_a_7 * (p_ImK22 + ((p_Eccentricity*p_Eccentricity) * ((3.0 * p_ImK10 / 4.0) + (p_ImK12 / 8.0) - (5.0 * p_ImK22) + (147.0 * p_ImK32 / 8.0)))); } /* @@ -798,9 +801,12 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.7) * * - * double CalculateDEccentricityTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDEccentricityTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * + * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_M2 Mass of companion star (Msol) @@ -809,13 +815,13 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in Eccentricity for binary (1/yr) */ - double CalculateDEccentricityTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDEccentricityTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_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; - return -(3.0 / 4.0) * (p_Eccentricity/p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2 / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * (-p_ImK22); + return -(3.0 / 4.0) * (p_Eccentricity/p_Omega) * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr * p_M2 / R1_AU / R1_AU / R1_AU) * R1_over_a_8 * ((3.0 * p_ImK10 / 2.0) - (p_ImK12 / 4.0) - p_ImK22 + (49.0 * p_ImK32 / 4.0)); } /* @@ -823,9 +829,12 @@ class BaseBinaryStar { * Zahn, 1977, Eq. (3.8) * * - * double CalculateDOmegaTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) + * double CalculateDOmegaTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) * + * @param [IN] p_ImK10 Imaginary (1,0) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK12 Imaginary (1,2) component of potential tidal love number of star (unitless) * @param [IN] p_ImK22 Imaginary (2,2) component of potential tidal love number of star (unitless) + * @param [IN] p_ImK32 Imaginary (3,2) component of potential tidal love number of star (unitless) * @param [IN] p_M1 Mass of star (Msol) * @param [IN] p_R1 Radius of star (Rsol) * @param [IN] p_I1 Moment of Inertia of star (Msol * AU^2) @@ -834,13 +843,13 @@ class BaseBinaryStar { * @param [IN] p_Eccentricity Eccentricity for binary * @return Change in Omega for star (1/yr/yr) */ - double CalculateDOmegaTidalDt(const double p_ImK22, const double p_M1, const double p_R1, const double p_I1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { + double CalculateDOmegaTidalDt(const double p_ImK10, const double p_ImK12, const double p_ImK22, const double p_ImK32, const double p_M1, const double p_R1, const double p_I1, const double p_M2, const double p_Omega, const double p_SemiMajorAxis, const double p_Eccentricity) { double R1_AU = p_R1 * RSOL_TO_AU; double R1_over_a = R1_AU / p_SemiMajorAxis; double R1_over_a_6 = R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a * R1_over_a; - return (3.0 / 2.0) * (1/p_I1) * (G_AU_Msol_yr * p_M2 * p_M2 / R1_AU) * R1_over_a_6 * p_ImK22; + return (3.0 / 2.0) * (1/p_I1) * (G_AU_Msol_yr * p_M2 * p_M2 / R1_AU) * R1_over_a_6 * (p_ImK22 + ((p_Eccentricity*p_Eccentricity) * ((p_ImK12 / 4.0) - (5.0 * p_ImK22) + (49.0 * p_ImK32 / 4.0)))); } }; diff --git a/src/BaseStar.h b/src/BaseStar.h index 1346b7571..f7e0bc3f8 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -179,7 +179,7 @@ class BaseStar { virtual void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { } // Default is NO-OP virtual void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables - virtual double CalculateImK22Tidal(const double p_Omega) { return 0.0; } // Default is 0.0 + virtual double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { return 0.0; } // Default is 0.0 void CalculateLambdas() { CalculateLambdas(m_Mass - m_CoreMass); } // Use class member variables void CalculateLambdas(const double p_EnvMass); diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index e60f9846e..f263339a1 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1131,16 +1131,18 @@ double GiantBranch::CalculateMomentOfInertia() const { } /* - * Calculate imaginary component of potential tidal Love number (Dynamical Tides, (2,2) mode ONLY) + * Calculate (n=2, l, m) imaginary component of potential tidal Love number (Dynamical Tides ONLY for now) * * Zahn, 1977, Eq. (5.5) , with the value of E_2 coming from Kushnir et al., 2017, by comparing Eq. (8) to Eq. (1) * - * double CalculateImK22Tidal(const double p_Omega) + * double CalculateImKlmTidal(const double p_Omega) * * @param [IN] p_Omega Orbital angular frequency (1/yr) - * @return Imaginary component of potential tidal love number (unitless, (2,2) mode only) + * @param [IN] l l mode + * @param [IN] m m mode + * @return Imaginary component of potential tidal love number (unitless) */ -double GiantBranch::CalculateImK22Tidal(const double p_Omega) { +double GiantBranch::CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { double beta2Dynamical = 1.0; double rhoFactorDynamcial = 0.1; double radiusAU = m_Radius * RSOL_TO_AU; @@ -1153,12 +1155,12 @@ double GiantBranch::CalculateImK22Tidal(const double p_Omega) { double E2Dynamical = (2.0 / 3.0) * coreRadius_over_radius_9 * mass_over_coreMass * std::cbrt(mass_over_coreMass) * beta2Dynamical * rhoFactorDynamcial; - double s22 = 2.0 * (p_Omega - Omega()) * std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass); - double s22_4_3 = s22 * std::cbrt(s22); - double s22_8_3 = s22_4_3 * s22_4_3; - double k22Dynamical = E2Dynamical * s22_8_3; + double slm = ((p_l*p_Omega) - (p_m*Omega())) * std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass); + double slm_4_3 = slm * std::cbrt(slm); + double slm_8_3 = slm_4_3 * slm_4_3; + double klmDynamical = E2Dynamical * slm_8_3; - return k22Dynamical; + return klmDynamical; } diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 0941f0a4a..ce50e5215 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -52,7 +52,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence { void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables static double CalculateHRateConstant_Static(const double p_Mass); - double CalculateImK22Tidal(const double p_Omega); + double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m); virtual double CalculateInitialSupernovaMass() const { return m_Mass; } // Use class member variables double CalculateLifetimeToHeIgnition(const double p_Mass, const double p_Tinf1_FGB, const double p_Tinf2_FGB); diff --git a/src/Star.h b/src/Star.h index 2b696e2fd..333adf533 100755 --- a/src/Star.h +++ b/src/Star.h @@ -162,7 +162,7 @@ class Star { double CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } double CalculateEddyTurnoverTimescale() { return m_Star->CalculateEddyTurnoverTimescale(); } - double CalculateImK22Tidal(const double p_Omega) { return m_Star->CalculateImK22Tidal(p_Omega); } + double CalculateImKlmTidal(const double p_Omega, const int p_l, const int p_m) { return m_Star->CalculateImKlmTidal(p_Omega, p_l, p_m); } void CalculateLambdas() { m_Star->CalculateLambdas(); } void CalculateLambdas(const double p_EnvMass) { m_Star->CalculateLambdas(p_EnvMass); }