Skip to content

Commit

Permalink
Implemented secular equations for tidal evolution following Zahn (1977)
Browse files Browse the repository at this point in the history
  • Loading branch information
@veome22
veome22 committed Apr 4, 2024
1 parent e7c35be commit 14ef044
Show file tree
Hide file tree
Showing 2 changed files with 135 additions and 12 deletions.
73 changes: 72 additions & 1 deletion src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2279,8 +2279,79 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) {

CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations

if (OPTIONS->EnableTides() && !m_Unbound) {
if (OPTIONS->EnableRealisticTides() && !m_Unbound) {
// Change binary semi-major axis 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),
m_Star1->Mass(),
m_Star1->Radius(),
m_Star2->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

double DSemiMajorAxis2Dt = CalculateDSemiMajorAxisTidalDt(m_Star2->CalculateImK22Tidal(m_Omega),
m_Star2->Mass(),
m_Star2->Radius(),
m_Star1->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

double DEccentricity1Dt = CalculateDEccentricityTidalDt(m_Star1->CalculateImK22Tidal(m_Omega),
m_Star1->Mass(),
m_Star1->Radius(),
m_Star2->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

double DEccentricity2Dt = CalculateDEccentricityTidalDt(m_Star2->CalculateImK22Tidal(m_Omega),
m_Star2->Mass(),
m_Star2->Radius(),
m_Star1->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

double DOmega1Dt = CalculateDOmegaTidalDt(m_Star1->CalculateImK22Tidal(m_Omega),
m_Star1->Mass(),
m_Star1->Radius(),
m_Star1->CalculateMomentOfInertiaAU(),
m_Star2->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

double DOmega2Dt = CalculateDOmegaTidalDt(m_Star2->CalculateImK22Tidal(m_Omega),
m_Star2->Mass(),
m_Star2->Radius(),
m_Star2->CalculateMomentOfInertiaAU(),
m_Star1->Mass(),
m_Omega,
m_SemiMajorAxis,
m_Eccentricity);

m_Star1->SetOmega(m_Star1->Omega() + (DOmega1Dt * p_Dt * MYR_TO_YEAR)); // synchronise star 1
m_Star2->SetOmega(m_Star2->Omega() + (DOmega2Dt * p_Dt * MYR_TO_YEAR)); // synchronise star 1

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_Omega = std::sqrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / m_SemiMajorAxis / m_SemiMajorAxis / m_SemiMajorAxis); // re-calculate orbital frequency

m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate total angular momentum

}

if (OPTIONS->EnableTides() && !m_Unbound) {
// find omega assuming synchronisation
// use current value of m_Omega as best guess for root
// if m_Omega == 0.0 (should only happen on the first timestep), calculate m_Omega here
Expand Down
74 changes: 63 additions & 11 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -767,25 +767,77 @@ class BaseBinaryStar {
}

/*
* Change in semi-major axis based on secular equations for tidal evolution given some tidal love number
* Change in semi-major axis based on secular equations for tidal evolution given the tidal love number
* Zahn, 1977, Eq. (3.6)
*
*
* double CalculateDSemiMajorAxisTidal(const double p_ImKlm, 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_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)
*
* @param [IN] p_ImKlm Imaginary component of potential tidal love number of star 1 (unitless)
* @param [IN] p_M1 Mass of star 1 (Msol)
* @param [IN] p_R1 Radius of star 1 (Rsol)
* @param [IN] p_M2 Mass of star 2 (Msol)
* @param [IN] p_Omega Orbital angular frequency for bianry (1/yr)
* @param [IN] p_ImK22 Imaginary (2,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)
* @param [IN] p_Omega Orbital angular frequency for binary (1/yr)
* @param [IN] p_SemiMajorAxis Semi-major axis for binary (AU)
* @param [IN] p_Eccentricity Eccentricity for binary
* @return Root found: will be -1.0 if no acceptable real root found
* @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 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;
}

/*
* 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)
*
* @param [IN] p_ImK22 Imaginary (2,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)
* @param [IN] p_Omega Orbital angular frequency for binary (1/yr)
* @param [IN] p_SemiMajorAxis Semi-major axis for binary (AU)
* @param [IN] p_Eccentricity Eccentricity for binary
* @return Change in Eccentricity for binary (1/yr)
*/
double CalculateDSemiMajorAxisTidal(const double p_ImKlm, 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 p_R1_AU = p_R1 * RSOL_TO_AU;
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 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);
}

/*
* 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)
*
* @param [IN] p_ImK22 Imaginary (2,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)
* @param [IN] p_M2 Mass of companion star (Msol)
* @param [IN] p_SemiMajorAxis Semi-major axis for binary (AU)
* @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 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 / p_Omega) * p_M2 * (1.0 + (p_M2 / p_M1)) * (G_AU_Msol_yr / p_R1_AU / p_R1_AU) * PPOW(p_R1_AU / p_SemiMajorAxis, 7) * p_ImKlm;
return (3.0 / 2.0) * (1/p_I1) * (G_AU_Msol_yr * p_M2 * p_M2 / R1_AU) * R1_over_a_6 * p_ImK22;
}

};
Expand Down

0 comments on commit 14ef044

Please sign in to comment.