diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index b87d93ea3..f01e21860 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -1,5 +1,5 @@ ##~!!~## COMPAS option values -##~!!~## File Created Tue Dec 9 15:41:35 2025 by COMPAS v03.27.02 +##~!!~## File Created Mon Feb 2 10:03:24 2026 by COMPAS v03.29.00 ##~!!~## ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the @@ -277,7 +277,7 @@ stringChoices: # --common-envelope-second-stage-gamma-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['ARBITRARY','KLENCKI_LINEAR','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS'] ### TIDES -# --tides-prescription: 'NONE' # Default: 'NONE' # Options: ['KAPIL2025','PERFECT','NONE'] +# --tides-prescription: 'NONE' # Default: 'NONE' # Options: ['ZAHN1977','KAPIL2026','PERFECT','NONE'] ### SUPERNOVAE, KICKS AND REMNANTS # --black-hole-kicks-mode: 'FALLBACK' # Default: 'FALLBACK' # Options: ['FALLBACK','ZERO','REDUCED','FULL'] diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index 88d23ca50..7229c67cc 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -1388,10 +1388,11 @@ A record is written to the System Snapshot logfile on the first timestep at whic **--tides-prescription** |br| Prescription for tidal evolution of the binary. |br| -Options: { NONE, PERFECT, KAPIL2025 } |br| +Options: { NONE, PERFECT, KAPIL2026, ZAHN1977 } |br| ``NONE`` disables tidal interactions. |br| ``PERFECT`` evolves the binary assuming instantaneous synchronization and circularization. |br| -``KAPIL2025`` uses the prescription from Kapil+ (2025). |br| +``KAPIL2026`` uses the prescription from Kapil+ (2026). |br| +``ZAHN1977`` is based on Zahn (1977) and Hurley+ (2002). |br| Default = NONE **--timestep-filename** |br| diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 898f91423..04887c1c5 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,10 @@ What's new Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``. +**03.29.00 January 30, 2026** +* Added new tidal prescription based on Zahn (1977) and Hurley et. al (2002), called ``ZAHN1977``. +* Renamed the ``KAPIL2025`` tides prescription to ``KAPIL2026`` to match paper date. + **03.28.00 January 28, 2026** * Updated the default values for the Mandel-Müller kick prescription to a magnitude of 630 km/s and a sigma of 0.45, as calibrated by Disberg et al. (2026) to the results of Disberg & Mandel (2025). diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 0f6a1cc78..83b551bf2 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1121,16 +1121,22 @@ double BaseBinaryStar::CalculateDOmegaTidalDt(const DBL_DBL_DBL_DBL p_ImKnm, con double radiusStar = p_Star->Radius(); double massCompanion = p_Star == m_Star1 ? m_Star2->Mass() : m_Star1->Mass(); + double e2 = m_Eccentricity * m_Eccentricity; + double e4 = e2 * e2; + double e6 = e4 * e2; + double ImK10, ImK12, ImK22, ImK32; std::tie(ImK10, ImK12, ImK22, ImK32) = p_ImKnm; 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_spin_term = (m_Eccentricity * m_Eccentricity) * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + double e2_spin_term = e2 * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + + double pseudo_sync_limit = OrbitalAngularVelocity() * (1.0 + (15.0 / 2.0) * e2 + (45.0 / 8.0) * e4 + (5.0 / 16.0) * e6) / ((1.0 + 3.0 * e2 + (3.0 / 8.0) * e4) * std::pow(1.0 - e2, 1.5)); // Hut 1981, Eq. (42) // 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_spin_term), 0.0) > 0)){e2_spin_term = 0.0;} + if ((utils::Compare(p_Star->Omega(), pseudo_sync_limit) > 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); } @@ -1152,19 +1158,25 @@ double BaseBinaryStar::CalculateDSemiMajorAxisTidalDt(const DBL_DBL_DBL_DBL p_Im double radiusStar = p_Star->Radius(); double massCompanion = p_Star == m_Star1 ? m_Star2->Mass() : m_Star1->Mass(); + double e2 = m_Eccentricity * m_Eccentricity; + double e4 = e2 * e2; + double e6 = e4 * e2; + double ImK10, ImK12, ImK22, ImK32; std::tie(ImK10, ImK12, ImK22, ImK32) = p_ImKnm; 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_sma_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 = e2 * ((3.0 * ImK10 / 4.0) + (ImK12 / 8.0) - (5.0 * ImK22) + (147.0 * ImK32 / 8.0)); + + double pseudo_sync_limit = OrbitalAngularVelocity() * (1.0 + (15.0 / 2.0) * e2 + (45.0 / 8.0) * e4 + (5.0 / 16.0) * e6) / ((1.0 + 3.0 * e2 + (3.0 / 8.0) * e4) * std::pow(1.0 - e2, 1.5)); // Hut 1981, Eq. (42) // 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;} + double e2_spin_term = e2 * ((ImK12 / 4.0) - (5.0 * ImK22) + (49.0 * ImK32 / 4.0)); + if ((utils::Compare(p_Star->Omega(), pseudo_sync_limit) > 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_sma_term); } @@ -2762,9 +2774,10 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { } break; - case TIDES_PRESCRIPTION::KAPIL2025: { // KAPIL2025 + case TIDES_PRESCRIPTION::KAPIL2026: { // KAPIL2026 - // Evolve binary semi-major axis, eccentricity, and spin of each star based on Kapil et al., 2025 + // Evolve binary semi-major axis, eccentricity, and spin of each star based on Kapil et al., 2026 + double jTotPrev = CalculateAngularMomentum(); DBL_DBL_DBL_DBL ImKnm1_tidal = m_Star1->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star2->Mass()); DBL_DBL_DBL_DBL ImKnm2_tidal = m_Star2->CalculateImKnmTidal(omega, m_SemiMajorAxis, m_Star1->Mass()); @@ -2784,12 +2797,23 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR))); fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR))); fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR))); - - m_Star1->SetOmega(m_Star1->Omega() + fraction_tidal_change * (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 1 spin - m_Star2->SetOmega(m_Star2->Omega() + fraction_tidal_change * (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 2 spin + m_SemiMajorAxis = m_SemiMajorAxis + fraction_tidal_change * ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve separation m_Eccentricity = m_Eccentricity + fraction_tidal_change * ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve eccentricity + + // correct stellar spins by an overall constant to enforce total angular momentum conservation while preserving spin ratio + // corrective_factor = (L_orb_i + I*omega1_i + I*omega2_i - L_orb_f) / (I1*omega1_f + I2*omega2_f) + // corrective_factor should be 1.0 if angular momentum is conserved properly by the secular equations + double angular_momentum_orbital = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + double omega1Proposed = m_Star1->Omega() + fraction_tidal_change * (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR); // evolve star 1 spin + double omega2Proposed = m_Star2->Omega() + fraction_tidal_change * (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR); // evolve star 2 spin + double omega_corrective_factor = (jTotPrev - angular_momentum_orbital) / (m_Star1->CalculateMomentOfInertiaAU() * omega1Proposed + m_Star2->CalculateMomentOfInertiaAU() * omega2Proposed); + omega_corrective_factor = std::isinf(omega_corrective_factor) || std::isnan(omega_corrective_factor) ? 1.0 : omega_corrective_factor; + + m_Star1->SetOmega(omega1Proposed * omega_corrective_factor); + m_Star2->SetOmega(omega2Proposed * omega_corrective_factor); + m_CircularizationTimescale = - m_Eccentricity / (DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * YEAR_TO_MYR; // Circularization timescale in Myr (for output files) m_CircularizationTimescale = (std::isnan(m_CircularizationTimescale) || std::isinf(m_CircularizationTimescale))? 0.0 : m_CircularizationTimescale; // check for NaN or Inf for circular binaries @@ -2849,6 +2873,133 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) { } } break; + case TIDES_PRESCRIPTION::ZAHN1977: { + // Evolve binary semi-major axis, eccentricity, and spin of each star based on Zahn (1977) + double envMass1, envMassMax1; + std::tie(envMass1, envMassMax1) = m_Star1->CalculateConvectiveEnvelopeMass(); + + double envMass2, envMassMax2; + std::tie(envMass2, envMassMax2) = m_Star2->CalculateConvectiveEnvelopeMass(); + + double pOrbital = _2_PI / omega; + double pSpin1 = _2_PI / m_Star1->Omega(); + double pSpin2 = _2_PI / m_Star2->Omega(); + double pTid1 = 1 / std::abs(1.0 / pOrbital - 1.0 / pSpin1); + double pTid2 = 1 / std::abs(1.0 / pOrbital - 1.0 / pSpin2); + + double tauConv1 = m_Star1->CalculateEddyTurnoverTimescale(); + double tauConv2 = m_Star2->CalculateEddyTurnoverTimescale(); + + double fConv1 = std::min(1.0, std::pow(pTid1 / (2.0 * tauConv1), 2)); + double fConv2 = std::min(1.0, std::pow(pTid2 / (2.0 * tauConv2), 2)); + double kOverT1 = (2.0 / 21.0) * (fConv1 / tauConv1) * (envMass1 / m_Star1->Mass()); + double kOverT2 = (2.0 / 21.0) * (fConv2 / tauConv2) * (envMass2 / m_Star2->Mass()); + + kOverT1 = std::isinf(kOverT1) ? 0.0 : kOverT1; // check for inf caused by zero tau_conv + kOverT2 = std::isinf(kOverT2) ? 0.0 : kOverT2; // check for inf caused by zero tau_conv + + double R1_AU = m_Star1->Radius() * RSOL_TO_AU; + double R2_AU = m_Star2->Radius() * RSOL_TO_AU; + + double omega1 = m_Star1->Omega(); + double omega2 = m_Star2->Omega(); + + double E2_Dynamical_Zahn1 = 1.592 * std::pow(10, -9) * std::pow(m_Star1->Mass(), 2.84); // From Hurley, et al. (2002), Eq. (43) + double E2_Dynamical_Zahn2 = 1.592 * std::pow(10, -9) * std::pow(m_Star2->Mass(), 2.84); // From Hurley, et al. (2002), Eq. (43) + + double w10 = 1.0 * omega; + + double w12_1 = (1.0 * omega - 2.0 * omega1); + double w12_2 = (1.0 * omega - 2.0 * omega2); + + double w22_1 = (2.0 * omega - 2.0 * omega1); + double w22_2 = (2.0 * omega - 2.0 * omega2); + + double w32_1 = (3.0 * omega - 2.0 * omega1); + double w32_2 = (3.0 * omega - 2.0 * omega2); + + // Compute equilibrium tidal Love numbers based on Zahn (1977), Eq. (4.1), (4.2) + double ImK22_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w22_1; + ImK22_Zahn_Equilibrium1 = std::isnan(ImK22_Zahn_Equilibrium1) ? 0.0 : ImK22_Zahn_Equilibrium1; // check for NaN + double ImK22_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w22_2; + ImK22_Zahn_Equilibrium2 = std::isnan(ImK22_Zahn_Equilibrium2) ? 0.0 : ImK22_Zahn_Equilibrium2; // check for NaN + + double ImK10_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w10; + ImK10_Zahn_Equilibrium1 = std::isnan(ImK10_Zahn_Equilibrium1) ? 0.0 : ImK10_Zahn_Equilibrium1; // check for NaN + double ImK10_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w10; + ImK10_Zahn_Equilibrium2 = std::isnan(ImK10_Zahn_Equilibrium2) ? 0.0 : ImK10_Zahn_Equilibrium2; // check for NaN + + double ImK12_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w12_1; + ImK12_Zahn_Equilibrium1 = std::isnan(ImK12_Zahn_Equilibrium1) ? 0.0 : ImK12_Zahn_Equilibrium1; // check for NaN + double ImK12_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w12_2; + ImK12_Zahn_Equilibrium2 = std::isnan(ImK12_Zahn_Equilibrium2) ? 0.0 : ImK12_Zahn_Equilibrium2; // check for NaN + + double ImK32_Zahn_Equilibrium1 = 2.0 * kOverT1 * (R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * w32_1; + ImK32_Zahn_Equilibrium1 = std::isnan(ImK32_Zahn_Equilibrium1) ? 0.0 : ImK32_Zahn_Equilibrium1; // check for NaN + double ImK32_Zahn_Equilibrium2 = 2.0 * kOverT2 * (R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * w32_2; + ImK32_Zahn_Equilibrium2 = std::isnan(ImK32_Zahn_Equilibrium2) ? 0.0 : ImK32_Zahn_Equilibrium2; // check for NaN + + // Compute dynamical tidal Love numbers based on Zahn (1977), Eq. (5.5) + double ImK22_Zahn_Dynamical1 = std::copysign(1.0, w22_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w22_1), 8.0/3.0); + ImK22_Zahn_Dynamical1 = std::isnan(ImK22_Zahn_Dynamical1) ? 0.0 : ImK22_Zahn_Dynamical1; // check for NaN + double ImK22_Zahn_Dynamical2 = std::copysign(1.0, w22_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w22_2), 8.0/3.0); + ImK22_Zahn_Dynamical2 = std::isnan(ImK22_Zahn_Dynamical2) ? 0.0 : ImK22_Zahn_Dynamical2; // check for NaN + + double ImK10_Zahn_Dynamical1 = std::copysign(1.0, w10) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w10), 8.0/3.0); + ImK10_Zahn_Dynamical1 = std::isnan(ImK10_Zahn_Dynamical1) ? 0.0 : ImK10_Zahn_Dynamical1; // check for NaN + double ImK10_Zahn_Dynamical2 =std::copysign(1.0, w10) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w10), 8.0/3.0); + ImK10_Zahn_Dynamical2 = std::isnan(ImK10_Zahn_Dynamical2) ? 0.0 : ImK10_Zahn_Dynamical2; // check for NaN + + double ImK12_Zahn_Dynamical1 = std::copysign(1.0, w12_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w12_1), 8.0/3.0); + ImK12_Zahn_Dynamical1 = std::isnan(ImK12_Zahn_Dynamical1) ? 0.0 : ImK12_Zahn_Dynamical1; // check for NaN + double ImK12_Zahn_Dynamical2 = std::copysign(1.0, w12_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w12_2), 8.0/3.0); + ImK12_Zahn_Dynamical2 = std::isnan(ImK12_Zahn_Dynamical2) ? 0.0 : ImK12_Zahn_Dynamical2; // check for NaN + + double ImK32_Zahn_Dynamical1 = std::copysign(1.0, w32_1) * E2_Dynamical_Zahn1 * std::pow(std::sqrt(R1_AU * R1_AU * R1_AU / (G_AU_Msol_yr * m_Star1->Mass())) * std::abs(w32_1), 8.0/3.0); + ImK32_Zahn_Dynamical1 = std::isnan(ImK32_Zahn_Dynamical1) ? 0.0 : ImK32_Zahn_Dynamical1; // check for NaN + double ImK32_Zahn_Dynamical2 = std::copysign(1.0, w32_2) * E2_Dynamical_Zahn2 * std::pow(std::sqrt(R2_AU * R2_AU * R2_AU / (G_AU_Msol_yr * m_Star2->Mass())) * std::abs(w32_2), 8.0/3.0); + ImK32_Zahn_Dynamical2 = std::isnan(ImK32_Zahn_Dynamical2) ? 0.0 : ImK32_Zahn_Dynamical2; // check for NaN + + // Combine equilibrium and dynamical tidal Love numbers + DBL_DBL_DBL_DBL ImKnm1_tidal = std::make_tuple(ImK10_Zahn_Equilibrium1+ImK10_Zahn_Dynamical1, ImK12_Zahn_Equilibrium1+ImK12_Zahn_Dynamical1, ImK22_Zahn_Equilibrium1 + ImK22_Zahn_Dynamical1, ImK32_Zahn_Equilibrium1 + ImK32_Zahn_Dynamical1); + DBL_DBL_DBL_DBL ImKnm2_tidal = std::make_tuple(ImK10_Zahn_Equilibrium2+ImK10_Zahn_Dynamical2, ImK12_Zahn_Equilibrium2+ImK12_Zahn_Dynamical2, ImK22_Zahn_Equilibrium2 + ImK22_Zahn_Dynamical2, ImK32_Zahn_Equilibrium2 + ImK32_Zahn_Dynamical2); + + // Perform secular evolution due to tides + double DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKnm1_tidal, m_Star1); // change in semi-major axis from star1 + double DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKnm2_tidal, m_Star2); // change in semi-major axis from star2 + + double DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(ImKnm1_tidal, m_Star1); // change in eccentricity from star1 + double DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(ImKnm2_tidal, m_Star2); // change in eccentricity from star2 + + double DOmega1Dt_tidal = CalculateDOmegaTidalDt(ImKnm1_tidal, m_Star1); // change in spin from star1 + double DOmega2Dt_tidal = CalculateDOmegaTidalDt(ImKnm2_tidal, m_Star2); // change in spin from star2 + + // limit change in stellar and orbital properties from tides to a maximum fraction of the current value + double fraction_tidal_change = 1.0; + fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR))); + fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * OrbitalAngularVelocity() / (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR))); + fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR))); + fraction_tidal_change = std::min(fraction_tidal_change, std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR))); + + m_Star1->SetOmega(m_Star1->Omega() + fraction_tidal_change * (DOmega1Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 1 spin + m_Star2->SetOmega(m_Star2->Omega() + fraction_tidal_change * (DOmega2Dt_tidal * p_Dt * MYR_TO_YEAR)); // evolve star 2 spin + m_SemiMajorAxis = m_SemiMajorAxis + fraction_tidal_change * ((DSemiMajorAxis1Dt_tidal + DSemiMajorAxis2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve separation + m_Eccentricity = m_Eccentricity + fraction_tidal_change * ((DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * p_Dt * MYR_TO_YEAR); // evolve eccentricity + + m_CircularizationTimescale = - m_Eccentricity / (DEccentricity1Dt_tidal + DEccentricity2Dt_tidal) * YEAR_TO_MYR; // Circularization timescale in Myr (for output files) + m_CircularizationTimescale = (std::isnan(m_CircularizationTimescale) || std::isinf(m_CircularizationTimescale))? 0.0 : m_CircularizationTimescale; // check for NaN or Inf for circular binaries + + m_SynchronizationTimescale1 = - (m_Star1->Omega() - omega) / DOmega1Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star1 in Myr (for output files) + m_SynchronizationTimescale1 = (std::isnan(m_SynchronizationTimescale1) || std::isinf(m_SynchronizationTimescale1))? 0.0 : m_SynchronizationTimescale1; // check for NaN or Inf for synchronized binaries + + m_SynchronizationTimescale2 = - (m_Star2->Omega() - omega) / DOmega2Dt_tidal * YEAR_TO_MYR; // Synchronization timescale for Star2 in Myr (for output files) + m_SynchronizationTimescale2 = (std::isnan(m_SynchronizationTimescale2) || std::isinf(m_SynchronizationTimescale2))? 0.0 : m_SynchronizationTimescale2; // check for NaN or Inf for synchronized binaries + + m_TotalAngularMomentum = CalculateAngularMomentum(); // re-calculate angular momenta + m_OrbitalAngularMomentum = CalculateOrbitalAngularMomentum(m_Star1->Mass(), m_Star2->Mass(), m_SemiMajorAxis, m_Eccentricity); + + } break; + default: // the only way this can happen is if someone added a TIDES_PRESCRIPTION // and it isn't accounted for in this code. We should not default here, with or without a warning. @@ -3065,7 +3216,7 @@ double BaseBinaryStar::ChooseTimestep(const double p_Factor) { dt = std::min(dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW); // yes - reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission } - if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::KAPIL2025) { // tides prescription = KAPIL2025 + if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::KAPIL2026) { // tides prescription = KAPIL2026 // yes - need to adjust dt double omega = OrbitalAngularVelocity(); @@ -3171,9 +3322,7 @@ void BaseBinaryStar::EvaluateBinary(const double p_Dt) { } CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations - ProcessTides(p_Dt); // process tides if required - // assign new values to "previous" values, for following timestep m_EccentricityPrev = m_Eccentricity; m_SemiMajorAxisPrev = m_SemiMajorAxis; diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 8220f1bee..629eb6703 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -3248,7 +3248,7 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do double radiusAU = m_Radius * RSOL_TO_AU; double coreRadiusAU = CalculateConvectiveCoreRadius() * RSOL_TO_AU; double convectiveEnvRadiusAU = CalculateRadialExtentConvectiveEnvelope() * RSOL_TO_AU; - double radiusIntershellAU = radiusAU - convectiveEnvRadiusAU; // Outer radial coordinate of radiative intershell + double radiusIntershellAU = radiusAU - convectiveEnvRadiusAU; // Outer radial coordinate of radiative intershell, or inner radial coordinate of convective envelope // There should be no Dynamical tides if the entire star is convective, i.e. if there are no convective-radiative boundaries. // If so, return 0.0 for all dynamical components of ImKnm. @@ -3289,34 +3289,34 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do double coreRadiusOverRadius_3 = coreRadiusOverRadius * coreRadiusOverRadius * coreRadiusOverRadius; double coreRadiusOverRadius_9 = coreRadiusOverRadius_3 * coreRadiusOverRadius_3 * coreRadiusOverRadius_3; double massOverCoreMass = m_Mass / coreMass; - double E2Dynamical = (2.0 / 3.0) * coreRadiusOverRadius_9 * massOverCoreMass * std::cbrt(massOverCoreMass) * beta2Dynamical * rhoFactorDynamcial; + double E2Core = (2.0 / 3.0) * coreRadiusOverRadius_9 * massOverCoreMass * std::cbrt(massOverCoreMass) * beta2Dynamical * rhoFactorDynamcial; // (l=2, n=1, m=0), Gravity Wave dissipation from core boundary double s10 = w10 * sqrtR3OverG_M; double s10_4_3 = s10 * std::cbrt(s10); double s10_8_3 = s10_4_3 * s10_4_3; - k10GravityCore = E2Dynamical * std::copysign(s10_8_3, w10); + k10GravityCore = E2Core * std::copysign(s10_8_3, w10); if (std::isnan(k10GravityCore)) k10GravityCore = 0.0; // (l=2, n=1, m=2), Gravity Wave dissipation from core boundary double s12 = w12 * sqrtR3OverG_M; double s12_4_3 = s12 * std::cbrt(s12); double s12_8_3 = s12_4_3 * s12_4_3; - k12GravityCore = E2Dynamical * std::copysign(s12_8_3, w12); + k12GravityCore = E2Core * std::copysign(s12_8_3, w12); if (std::isnan(k12GravityCore)) k12GravityCore = 0.0; // (l=2, n=2, m=2), Gravity Wave dissipation from core boundary double s22 = w22 * sqrtR3OverG_M; double s22_4_3 = s22 * std::cbrt(s22); double s22_8_3 = s22_4_3 * s22_4_3; - k22GravityCore = E2Dynamical * std::copysign(s22_8_3, w22); + k22GravityCore = E2Core * std::copysign(s22_8_3, w22); if (std::isnan(k22GravityCore)) k22GravityCore = 0.0; // (l=2, n=3, m=2), Gravity Wave dissipation from core boundary double s32 = w32 * sqrtR3OverG_M; double s32_4_3 = s32 * std::cbrt(s32); double s32_8_3 = s32_4_3 * s32_4_3; - k32GravityCore = E2Dynamical * std::copysign(s32_8_3, w32); + k32GravityCore = E2Core * std::copysign(s32_8_3, w32); if (std::isnan(k32GravityCore)) k32GravityCore = 0.0; } @@ -3328,8 +3328,8 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do if ((utils::Compare(convectiveEnvRadiusAU / radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) > 0) || (utils::Compare(envMass / m_Mass, TIDES_MINIMUM_FRACTIONAL_EXTENT) > 0)) { constexpr double dynPrefactor = 3.207452512782476; // 3^(11/3) * Gamma(1/3)^2 / 40 PI - constexpr double m_l_factor_22 = 0.183440402716368; // m * (l(l+1))^{-4/3} - double cbrtdNdlnr = std::cbrt(G_AU_Msol_yr * radIntershellMass / radiusIntershellAU / (radiusAU - radiusIntershellAU) / (radiusAU - radiusIntershellAU)); + constexpr double m_l_factor_22 = 0.091720201358184; // (l(l+1))^{-4/3}, assuming l=2 + double cbrtdNdlnr = std::cbrt(G_AU_Msol_yr * radIntershellMass / radiusIntershellAU / radiusIntershellAU / (radiusAU - radiusIntershellAU)); double alpha = radiusIntershellAU / radiusAU; double oneMinusAlpha = 1.0 - alpha; @@ -3349,27 +3349,33 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmDynamical(const double p_Omega, const do double alpha_2_3Minus_1 = (alpha * 2.0 / 3.0) - 1.0; // Assume GW dissipation from the envelope boundary only acts if the radiative zone extends to the core, i.e. if there is no convective core. - if (utils::Compare(coreRadiusAU/radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) < 0) { + if (utils::Compare(coreRadiusAU/radiusAU, TIDES_MINIMUM_FRACTIONAL_EXTENT) < 0) { double Epsilon = alpha_11 * envMass / m_Mass * oneMinusGamma_2 * alpha_2_3Minus_1 * alpha_2_3Minus_1 / beta_2 / oneMinusAlpha_3 / oneMinusAlpha_2; + double R3_Epsilon_dNdlnr_factor = R3OverG_M * Epsilon / cbrtdNdlnr; + double E2Envelope = dynPrefactor * m_l_factor_22 * R3_Epsilon_dNdlnr_factor; - // (l=2, n=1, m=0), Gravity Wave dissipation from envelope boundary is always 0.0 since m * (l(l+1))^{-4/3} = 0 + // (l=2, n=1, m=0), Gravity Wave dissipation from envelope boundary + double w10_4_3 = w10 * std::cbrt(w10); + double w10_8_3 = w10_4_3 * w10_4_3; + k10GravityEnv = E2Envelope * std::copysign(w10_8_3, w10); + if (std::isnan(k10GravityEnv)) k10GravityEnv = 0.0; // (l=2, n=1, m=2), Gravity Wave dissipation from envelope boundary double w12_4_3 = w12 * std::cbrt(w12); double w12_8_3 = w12_4_3 * w12_4_3; - k12GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w12_8_3, w12) * R3OverG_M * Epsilon / cbrtdNdlnr; + k12GravityEnv = E2Envelope * std::copysign(w12_8_3, w12); if (std::isnan(k12GravityEnv)) k12GravityEnv = 0.0; // (l=2, n=2, m=2), Gravity Wave dissipation from envelope boundary double w22_4_3 = w22 * std::cbrt(w22); double w22_8_3 = w22_4_3 * w22_4_3; - k22GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w22_8_3, w22)* R3OverG_M * Epsilon / cbrtdNdlnr; + k22GravityEnv = E2Envelope * std::copysign(w22_8_3, w22); if (std::isnan(k22GravityEnv)) k22GravityEnv = 0.0; // (l=2, n=3, m=2), Gravity Wave dissipation from envelope boundary double w32_4_3 = w32 * std::cbrt(w32); double w32_8_3 = w32_4_3 * w32_4_3; - k32GravityEnv = dynPrefactor * m_l_factor_22 * std::copysign(w32_8_3, w32) * R3OverG_M * Epsilon / cbrtdNdlnr; + k32GravityEnv = E2Envelope * std::copysign(w32_8_3, w32); if (std::isnan(k32GravityEnv)) k32GravityEnv = 0.0; } @@ -3436,7 +3442,7 @@ DBL_DBL_DBL_DBL BaseStar::CalculateImKnmEquilibrium(const double p_Omega, const double twoOmegaSpin = omegaSpin + omegaSpin; double rhoConv = envMass / (4.0 * M_PI * (rOut_3 - rIn_3) / 3.0); - double lConv = rEnvAU; // set length scale to height of convective envelope + double lConv = rEnvAU / 2.0; // set length scale to height of convective envelope double tConv = CalculateEddyTurnoverTimescale(); double vConv = lConv / tConv; double omegaConv = 1.0 / tConv; // absent factor of 2*PI, following Barker (2020) diff --git a/src/Options.h b/src/Options.h index 39a19fd01..b22eafe02 100755 --- a/src/Options.h +++ b/src/Options.h @@ -241,6 +241,7 @@ class Options { { "pulsational-pair-instability-prescription", "COMPAS", "WOOSLEY", false, "20250208" }, { "pulsar-birth-spin-period-distribution", "ZERO", "NOSPIN", false, "20250303" }, { "tides-prescription", "KAPIL2024", "KAPIL2025", false, "20250525" }, + { "tides-prescription", "KAPIL2025", "KAPIL2026", false, "20260130" }, { "mass-loss-prescription", "MERRITT2024", "MERRITT2025", false, "20250717" }, { "use-mass-loss", "TRUE", "MERRITT2025", true, "20250809" }, { "use-mass-loss", "ON", "MERRITT2025", true, "20250809" }, diff --git a/src/changelog.h b/src/changelog.h index 67ffd087a..0ff319a51 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1688,7 +1688,12 @@ // - Fix issue #1446: Theta and phi variables are flipped when assigning kicks, potentially giving unintended kick distributions // 03.28.00 PD - January 28, 2026 - Enhancement: // - Updated default MullerMandel kick parameters in constants.h to the values from Disberg+2026, previous values were from Kapil+2023 -// +// 03.29.00 VK - January 30, 2026 - Enhancements: +// - Added new tidal prescription based on Zahn (1977) and Hurley et. al (2002), called '--tides-prescription ZAHN1977'. +// - Renamed KAPIL2025 tides prescription to KAPIL2026 to match publication. +// - Updated the spin limit in 'BaseBinaryStar::CalculateDOmegaTidalDt()' to allow pseudo-synchronization based on Hut (1981), which affects the maximum spin with KAPIL2026 and ZAHN1977 options. +// - Added a limit to rotation change per time step in KAPIL2026 to ensure angular momentum conservation. +// - Updated dynamical tides equations in KAPIL2026 to match paper. // // Version string format is MM.mm.rr, where // @@ -1699,7 +1704,7 @@ // if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made // if mm is incremented, set rr to 00, even if defect repairs were also made -const std::string VERSION_STRING = "03.28.00"; +const std::string VERSION_STRING = "03.29.00"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index 53880af4c..ddfb8ec3e 100755 --- a/src/constants.h +++ b/src/constants.h @@ -298,7 +298,7 @@ constexpr int TIDES_OMEGA_MAX_TRIES = 30; constexpr int TIDES_OMEGA_MAX_ITERATIONS = 50; // Maximum number of root finder iterations in BaseBinaryStar::OmegaAfterCircularisation() constexpr double TIDES_OMEGA_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in BaseBinaryStar::OmegaAfterCircularisation() (added to 1.0) constexpr double TIDES_MINIMUM_FRACTIONAL_EXTENT = 1.0E-4; // Minimum fractional radius or mass of the stellar core or envelope, above which a given tidal dissipation mechanism is considered applicable -constexpr double TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC = 0.01; // Maximum allowed change in orbital and spin properties due to KAPIL2025 tides in a single timestep - 1% expressed as a fraction +constexpr double TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC = 0.01; // Maximum allowed change in orbital and spin properties due to KAPIL2026 tides in a single timestep - 1% expressed as a fraction constexpr double TIDES_MINIMUM_FRACTIONAL_NUCLEAR_TIME = 0.001; // Minimum allowed timestep from tidal processes, as a fraction of the nuclear minimum time scale constexpr double FARMER_PPISN_UPP_LIM_LIN_REGIME = 38.0; // Maximum CO core mass to result in the linear remnant mass regime of the FARMER PPISN prescription diff --git a/src/typedefs.h b/src/typedefs.h index 04a2f9690..7d16ddbe9 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -937,11 +937,12 @@ const COMPASUnorderedMap STELLAR_POPULATION_LAB }; // tides prescriptions -enum class TIDES_PRESCRIPTION: int { NONE, PERFECT, KAPIL2025 }; +enum class TIDES_PRESCRIPTION: int { NONE, PERFECT, KAPIL2026, ZAHN1977 }; const COMPASUnorderedMap TIDES_PRESCRIPTION_LABEL = { { TIDES_PRESCRIPTION::NONE, "NONE" }, { TIDES_PRESCRIPTION::PERFECT, "PERFECT" }, - { TIDES_PRESCRIPTION::KAPIL2025, "KAPIL2025" } + { TIDES_PRESCRIPTION::KAPIL2026, "KAPIL2026" }, + { TIDES_PRESCRIPTION::ZAHN1977, "ZAHN1977" } }; // symbolic names for timescales