Major changes to tidal prescriptions

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']

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|

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).

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)

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" },

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__

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

src/typedefs.h

@@ -937,11 +937,12 @@ const COMPASUnorderedMap<STELLAR_POPULATION, std::string> 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, std::string> 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