Skip to content

Commit

Permalink
Generalized tidal love number and resulting secular equations to othe…
Browse files Browse the repository at this point in the history 
…r (l,m) modes
  • Loading branch information
@veome22
veome22 committed Apr 8, 2024
1 parent b62e20e commit fd038ce
Show file tree
Hide file tree
Showing 6 changed files with 46 additions and 29 deletions.
22 changes: 14 additions & 8 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2335,47 +2335,52 @@ 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) {
m_Omega = OrbitalAngularVelocity();
}


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(),
m_Omega,
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(),
m_Omega,
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(),
m_Omega,
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(),
m_Omega,
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(),
Expand All @@ -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(),
Expand All @@ -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

Expand Down
27 changes: 18 additions & 9 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -785,22 +788,25 @@ 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))));
}

/*
* Change in eccentricity based on secular equations for tidal evolution given the tidal love number
* 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)
Expand All @@ -809,23 +815,26 @@ 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));
}

/*
* Change in spin based on secular equations for tidal evolution given the tidal love number
* 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)
Expand All @@ -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))));
}

};
Expand Down
2 changes: 1 addition & 1 deletion src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
20 changes: 11 additions & 9 deletions src/GiantBranch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
}


Expand Down
2 changes: 1 addition & 1 deletion src/GiantBranch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/Star.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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); }

Expand Down

0 comments on commit fd038ce

Please sign in to comment.