Add gravitational radiation emission to binary evolution.

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -14,6 +14,7 @@ booleanChoices:
 #    --enable-warnings: False                                              # Default: False                                                                                                                                                              # option to enable/disable warning messages
 #    --errors-to-file: False                                               # Default: False
 #    --evolve-unbound-systems: True                                        # Default: True
+#    --emit-gravitational-radiation: True                                  # Default: True
 #    --population-data-printing: False                                     # Default: False
 #    --print-bool-as-string: False                                         # Default: False
 #    --quiet: False                                                        # Default: False

online-docs/pages/User guide/Program options/program-options-list-defaults.rst

@@ -349,6 +349,10 @@ Default = 0.0
 Multiplication factor for Eddington accretion for NS & BH (i.e. > 1 is super-eddington and 0 is no accretion). |br|
 Default = 1.0
 
+**--emit-gravitational-radiation**  |br|
+Emit gravitational radiation at each timestep of binary evolution according to Peters 1964. |br|
+Default = TRUE
+
 **--enable-warnings** |br|
 Display warning messages to stdout. |br|
 Default = FALSE
@@ -1370,7 +1374,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 
 **Administrative**
 
---mode, --number-of-systems, --evolve-double-white-dwarfs, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction, --maximum-evolution-time, --maximum-number-timestep-iterations,
+--mode, --number-of-systems, --emit-gravitational-radiation, --evolve-double-white-dwarfs, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction, --maximum-evolution-time, --maximum-number-timestep-iterations,
 --radial-change-fraction, --random-seed, --timestep-multiplier, --timestep-filename
 
 --grid, --grid-start-line, --grid-lines-to-process

online-docs/pages/whats-new.rst

@@ -6,6 +6,11 @@ Following is a brief list of important updates to the COMPAS code.  A complete r
 
 **LATEST RELEASE** |br|
 
+**02.47.00 May 17, 2024**
+
+* New option to emit gravitational radiation at each timestep of binary evolution: ``--emit-gravitational-radiation``. The effects of radiation are approximated by the change in semimajor axis and eccentricity from Peters 1946 equations 5.6 and 5.7.
+* If ``--emit-gravitational-radiation``, we update timesteps dynamically as a function of :math:`da/dt` if gravitational radiation is large enough.
+
 **02.46.00 May 13, 2024**
 
 * added options ``--radial-change-fraction`` and ``--mass-change-fraction``, as approximate desired fractional changes in stellar radius and mass on phase when setting SSE and BSE timesteps

src/BaseBinaryStar.cpp

@@ -2370,6 +2370,69 @@ void BaseBinaryStar::ResolveMassChanges() {
 }
 
 
+/*
+ * Update the GW effects on the binary
+ *
+ * Use Peters 1964's approximations of the effects of emitting a GW by:
+ * - Updating the binary's semi-major axis (eq 5.6)
+ * - Updating the binary's eccentricity (eq 5.7)
+ * 
+ * void CalculateGravitationalRadiation()
+ * 
+ * @return                  The change in semi major axis and eccentricity per time due to GW emission
+*/
+std::tuple<double, double> BaseBinaryStar::CalculateGravitationalRadiation() {
+
+    // Useful values
+    double eccentricitySquared = m_Eccentricity * m_Eccentricity;
+    double oneMinusESq = 1.0 - eccentricitySquared;
+    double oneMinusESq_5 = oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq;
+    double G_AU_Msol_yr_3 = G_AU_Msol_yr * G_AU_Msol_yr * G_AU_Msol_yr;
+    double C_AU_Yr_5 = C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr;
+    double m_SemiMajorAxis_3 = m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis;
+    double massAndGAndCTerm = G_AU_Msol_yr_3 * m_Star1->Mass() * m_Star2->Mass() * (m_Star1->Mass() + m_Star2->Mass()) / C_AU_Yr_5;  // G^3 * m1 * m2(m1 + m2) / c^5 in units of Msol, AU and yr
+
+    // Approximate rate of change in semimajor axis
+    double numeratorA = -64.0 * massAndGAndCTerm;
+    double denominatorA = 5.0 * m_SemiMajorAxis_3 * std::sqrt(oneMinusESq_5 * oneMinusESq * oneMinusESq);
+    double DaDtGW = (numeratorA / denominatorA) * (1.0 + (73.0 / 24.0) * eccentricitySquared + (37.0 / 96.0) * eccentricitySquared * eccentricitySquared) * MYR_TO_YEAR;  // units of AU Myr^-1
+
+    // Approximate rate of change in eccentricity
+    double numeratorE = -304.0 * m_Eccentricity * massAndGAndCTerm;
+    double denominatorE = 15.0 * m_SemiMajorAxis_3 * m_SemiMajorAxis * std::sqrt(oneMinusESq_5);
+    double DeDtGW = (numeratorE / denominatorE) * (1.0 + (121.0 / 304.0) * eccentricitySquared) * YEAR_TO_MYR;  // units of Myr^-1
+
+    return std::make_tuple(DaDtGW, DeDtGW);
+}
+
+
+/*
+ * Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation().
+ * 
+ * This function updates the sem-major axis, eccentricity, and previous eccentricity values
+ * (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev) as a result of emitting the GW.
+ * 
+ * void EmitGravitationalRadiation(const double p_Dt)
+ *
+ * @param   [IN]    p_Dt                        timestep in Myr
+ * @param   [IN]    DaDtGW                      change in semimajor axis per time in AU/Myr due to GW radiation
+ * @param   [IN]    DeDtGW                      change in eccentricity per time in Myr^-1 due to GW radiation
+*/
+void BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW) {
+
+    // Update semimajor axis
+    double aNew = m_SemiMajorAxis + (DaDtGW * p_Dt);
+    m_SemiMajorAxis = aNew < 0.0 ? 1E-20 : (aNew);  // if <0, set to arbitrarily small number
+
+    // Update the eccentricity
+    m_Eccentricity += DeDtGW * p_Dt;
+
+    // Save values as previous timestep
+    m_SemiMajorAxisPrev = m_SemiMajorAxis;
+    m_EccentricityPrev = m_Eccentricity;
+}
+
+
 /*
  * Evaluate the binary system
  *
@@ -2674,6 +2737,9 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
             dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                                  // quantised
         }
 
+        double DaDtGW;
+        double DeDtGW;
+
         unsigned long int stepNum = 1;                                                                                                      // initialise step number
 
         while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                             // perform binary evolution - iterate over timesteps until told to stop
@@ -2708,6 +2774,10 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
 
                 if (OPTIONS->RLOFPrinting()) StashRLOFProperties(MASS_TRANSFER_TIMING::PRE_MT);                                             // stash properties immediately pre-Mass Transfer 
 
+                if (OPTIONS->EmitGravitationalRadiation()) {                                                                                // emit gravitational wave if desired
+                    EmitGravitationalWave(dt, DaDtGW, DeDtGW);                                                                              // this may update m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev
+                }
+
                 EvaluateBinary(dt);                                                                                                         // evaluate the binary at this timestep
 
                 (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_BINARY_TIMESTEP);                                            // print (log) detailed output
@@ -2793,6 +2863,11 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
 
             (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED);                                                  // print (log) detailed output: this is after all changes made in the timestep
 
+            if (evolutionStatus == EVOLUTION_STATUS::STELLAR_MERGER || evolutionStatus == EVOLUTION_STATUS::STARS_TOUCHING) {               // are the stars touching/merging?
+                m_MassTransferTrackerHistory = MT_TRACKING::MERGER;                                                                         // set flag for merger
+                (void)PrintCommonEnvelope();                                                                                                // print (log) common envelope details
+            }
+
             if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP;                          // number of timesteps for evolution exceeds maximum
             else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) {
                 evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED;                                                                    // using user-provided timesteps and all consumed
@@ -2814,6 +2889,13 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                 else {                                                                                                                      // not using user-provided timesteps
                     dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier();              // calculate new timestep
                     dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                          // quantised
+
+                    if (OPTIONS->EmitGravitationalRadiation()) {                                                                            // if our binary is emitting GWs
+                        std::tie(DaDtGW, DeDtGW) = CalculateGravitationalRadiation();                                                       // update the GW radiation effects da/dt and de/dt
+                        if (!usingProvidedTimesteps) {
+                            dt = std::min(dt, std::round(-1E-2 * m_SemiMajorAxis / DaDtGW / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM);          // make dt smaller if strong GW effects
+                        }
+                    }
                 }
 
                 if ((m_Star1->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) || m_Star2->IsOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) || dt < NUCLEAR_MINIMUM_TIMESTEP) {

src/BaseBinaryStar.h

@@ -416,6 +416,9 @@ class BaseBinaryStar {
 
     double  CalculateAngularMomentum() const                                    { return CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU()); }
 
+    std::tuple<double, double>    CalculateGravitationalRadiation();
+    void    EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW);
+
     void    CalculateEnergyAndAngularMomentum();
 
     double CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL 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);

src/Options.cpp

@@ -161,7 +161,8 @@ void Options::OptionValues::Initialise() {
     m_EvolveDoubleWhiteDwarfs                                       = false;
     m_EvolvePulsars                                                 = false;
 	m_EvolveUnboundSystems                                          = true;
-
+    m_EmitGravitationalRadiation                                    = true;
+
     m_NatalKickForPPISN                                             = false;
 
     m_DetailedOutput                                                = false;
@@ -807,6 +808,11 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt
             po::value<bool>(&p_Options->m_EvolveUnboundSystems)->default_value(p_Options->m_EvolveUnboundSystems)->implicit_value(true),                                                          
             ("Continue evolving stars even if the binary is disrupted (default = " + std::string(p_Options->m_EvolveUnboundSystems ? "TRUE" : "FALSE") + ")").c_str()
         )
+        (
+            "emit-gravitational-radiation",                                      
+            po::value<bool>(&p_Options->m_EmitGravitationalRadiation)->default_value(p_Options->m_EmitGravitationalRadiation)->implicit_value(true),                                                          
+            ("Emit gravitational radiation at each timestep of binary evolution (default = " + std::string(p_Options->m_EmitGravitationalRadiation ? "TRUE" : "FALSE") + ")").c_str()
+        )
         (
             "expel-convective-envelope-above-luminosity-threshold",
             po::value<bool>(&p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->default_value(p_Options->m_ExpelConvectiveEnvelopeAboveLuminosityThreshold)->implicit_value(true),

src/Options.h

@@ -336,6 +336,7 @@ class Options {
         "evolve-double-white-dwarfs",
         "evolve-pulsars",
         "evolve-unbound-systems",
+        "emit-gravitational-radiation",
 
         "initial-mass-1",
         "initial-mass-2",
@@ -450,6 +451,7 @@ class Options {
         "evolve-double-white-dwarfs"
         "evolve-pulsars",
         "evolve-unbound-systems",
+        "emit-gravitational-radiation",
 
         "fryer-supernova-engine",
 
@@ -674,6 +676,7 @@ class Options {
             bool                                                m_EvolvePulsars;                                                // Whether to evolve pulsars or not
             bool                                                m_NatalKickForPPISN;                                            // Flag if PPISN remnant should receive a non-zero natal kick
 	        bool                                                m_EvolveUnboundSystems;							                // Option to chose if unbound systems are evolved until death or the evolution stops after the system is unbound during a SN.
+            bool                                                m_EmitGravitationalRadiation;                                   // Option to emit gravitational radiation for each timestep of binary evolution
 
             bool                                                m_DetailedOutput;                                               // Print detailed output details to file (default = false)
             bool                                                m_PopulationDataPrinting;                                       // Print certain data for small populations, but not for larger one
@@ -1265,6 +1268,7 @@ class Options {
     bool                                        EvolveDoubleWhiteDwarfs() const                                         { return OPT_VALUE("evolve-double-white-dwarfs", m_EvolveDoubleWhiteDwarfs, true); }
     bool                                        EvolvePulsars() const                                                   { return OPT_VALUE("evolve-pulsars", m_EvolvePulsars, true); }
     bool                                        EvolveUnboundSystems() const                                            { return OPT_VALUE("evolve-unbound-systems", m_EvolveUnboundSystems, true); }
+    bool                                        EmitGravitationalRadiation() const                                      { return OPT_VALUE("emit-gravitational-radiation", m_EmitGravitationalRadiation, true); }
     bool                                        ExpelConvectiveEnvelopeAboveLuminosityThreshold() const                 { return OPT_VALUE("expel-convective-envelope-above-luminosity-threshold", m_ExpelConvectiveEnvelopeAboveLuminosityThreshold, true); }
 
     bool                                        FixedRandomSeedCmdLine() const                                          { return m_CmdLine.optionValues.m_FixedRandomSeed; }

src/changelog.h

@@ -1181,8 +1181,11 @@
 //                                      - fix for issue #744 - GB parameters `p` and `q` calculated differently for naked helium stars (see issue for details)
 //                                      - changed name of `ResolveEnvelopeLoss()` parameter `p_NoCheck` to `p_Force` (it is required, and now we understand why... see issue #873)
 //                                      - some code cleanup
+// 02.47.00    APB - May 17, 2024  = Enhancement:
+//                                      - Implemented gravitational radiation at each timestep of binary evolution. Available with new '--emit-gravitational-radiation' option.
+//                                          - If --emit-gravitational-radiation, timestep dynamically as function of da/dt.
 
 
-const std::string VERSION_STRING = "02.46.05";
+const std::string VERSION_STRING = "02.47.00";
 
 # endif // __changelog_h__

src/constants.h

@@ -273,6 +273,7 @@ constexpr double LSOLW                                  = 3.844E26;
 constexpr double AU                                     = 149597870700.0;                                           // 1 AU (Astronomical Unit) in metres
 constexpr double KM                                     = 1000.0;                                                   // 1 km (Kilometre) in metres
 constexpr double C                                      = 3.0E8;                                                    // Speed of light in m s^-1
+constexpr double C_AU_yr                                = C * 0.001 * KM_TO_AU * SECONDS_IN_YEAR;                   // Speed of light in AU yr^-1
 
 constexpr double MU_0                                   = 4.0 * M_PI * 1.0E-7;                                      // Vacuum permeability in m kg s-2 A-2
 

src/yaml.h

@@ -80,6 +80,7 @@ namespace yaml {
             "    --enable-warnings                                               # option to enable/disable warning messages",
             "    --errors-to-file",
             "    --evolve-unbound-systems",
+            "    --emit-gravitational-radiation",
             "    --population-data-printing",
             "    --print-bool-as-string",
             "    --quiet",