Skip to content

Commit

Permalink
Add gravitational radiation emission to binary evolution.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Adam-Boesky
Adam-Boesky committed Aug 23, 2024
1 parent 21b1eef commit 6a9fa62
Show file tree
Hide file tree
Showing 10 changed files with 158 additions and 45 deletions.
1 change: 1 addition & 0 deletions compas_python_utils/preprocessing/compasConfigDefault.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ booleanChoices:
# --errors-to-file: False # Default: False
# --evolve-main-sequence-mergers: False # Default: False
# --evolve-unbound-systems: True # Default: True
# --emit-gravitational-radiation: False # Default: False
# --population-data-printing: False # Default: False
# --print-bool-as-string: False # Default: False
# --quiet: False # Default: False
Expand Down
6 changes: 5 additions & 1 deletion online-docs/pages/User guide/Program options/program-options-list-defaults.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,6 +379,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 = FALSE

**--enable-warnings** |br|
Display warning messages to stdout. |br|
Default = FALSE
Expand Down Expand Up @@ -1440,7 +1444,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-main-sequence-mergers, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction,
--mode, --number-of-systems, --emit-gravitational-radiation, --evolve-double-white-dwarfs, --evolve-main-sequence-mergers, --evolve-pulsars, --evolve-unbound-systems, --mass-change-fraction,
--maximum-evolution-time, --maximum-number-timestep-iterations,
--radial-change-fraction, --random-seed, --timestep-multiplier, --timestep-filename

Expand Down
4 changes: 4 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@ Following is a brief list of important updates to the COMPAS code. A complete r

**LATEST RELEASE** |br|

**03.01.00 Aug 20, 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. Reduce timestep if required to keep orbital separation change per step due to GW radiation within ~ 1%.

**03.00.00 Jul 26, 2024**

This is a major release of COMPAS. There are some significant changes in COMPAS operation and functionality in this release. The major change, and the impetus for
Expand Down
152 changes: 117 additions & 35 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2490,6 +2490,93 @@ void BaseBinaryStar::ProcessTides(const double p_Dt) {
}


/*
* Calculate and emit gravitational radiation.
*
* This function uses Peters 1964 to approximate the effects of GW emission with two steps:
* - Calculate the change in semi-major axis (m_SemiMajorAxis) per time given by eq 5.6.
* - Calculate the change in eccentricity (m_Eccentricity) per time given by eq 5.7.
*
* m_DaDtGW and m_DeDtGW are updated so that they can be used to calculate the timestep dynamically.
*
* void CalculateGravitationalRadiation()
*/
void 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);
m_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);
m_DeDtGW = (numeratorE / denominatorE) * (1.0 + (121.0 / 304.0) * eccentricitySquared) * YEAR_TO_MYR; // units of Myr^-1
}


/*
* Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation().
*
* This function updates the semi-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
*/
void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {

// Update semimajor axis
double aNew = m_SemiMajorAxis + (m_DaDtGW * p_Dt);
m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number

// Update the eccentricity
m_Eccentricity += m_DeDtGW * p_Dt;

// Save values as previous timestep
m_SemiMajorAxisPrev = m_SemiMajorAxis;
m_EccentricityPrev = m_Eccentricity;
}


/*
* Choose a timestep based on the parameters of the binary.
*
* This function will return the minimum of (i) a timestep based on the
* orbital timescale of the binary and (ii) (if configured to emit GWs)
* a timestep based on the magnitude of gravitational radiation.
*
* double ChooseTimestep(double p_Dt)
*
* @param [IN] p_Dt previous timestep in Myr
* @return new timestep in Myr
*/
double BaseBinaryStar::ChooseTimestep(double p_Dt) {

p_Dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()); // calculate new timestep

if (OPTIONS->EmitGravitationalRadiation()) { // emitting GWs?
p_Dt = std::min(p_Dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW); // reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
}

p_Dt *= OPTIONS->TimestepMultiplier();

p_Dt = std::max(std::round(p_Dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum

return p_Dt;
}


/*
* Evaluate the binary system
*
Expand Down Expand Up @@ -2692,23 +2779,40 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
// evolve the current binary up to the maximum evolution time (and number of steps)

double dt; // timestep
if (usingProvidedTimesteps) { // user-provided timesteps?
// get new timestep
// - don't quantise
// - don't apply timestep multiplier
// (we assume user wants the timesteps in the file)
dt = timesteps[0];
}
else { // no - not using user-provided timesteps
dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier() / 1000.0; // calculate timestep - make first step small
dt = std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP); // quantised
}

unsigned long int stepNum = 1; // initialise step number
unsigned long int stepNum = 0; // initialise step number

while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // perform binary evolution - iterate over timesteps until told to stop

stepNum++; // increment stepNum

// if user selects to emit GWs, calculate the effects of radiation
// - note that this is placed before ChooseTimestep() is called because
// the timestep is a function of graviational radiation
if (OPTIONS->EmitGravitationalRadiation()) {
CalculateGravitationalRadiation();
}

if (stepNum > 1) { // after the first timestep, set previous timestep
m_Star2->UpdatePreviousTimestepDuration();
m_Star1->UpdatePreviousTimestepDuration();
}
if (usingProvidedTimesteps) { // user-provided timesteps?
// select a timestep
// - don't quantise
// - don't apply timestep multiplier
// (we assume user wants the timesteps in the file)
dt = timesteps[stepNum - 1];
}
else { // no - not using user-provided timesteps
dt = ChooseTimestep(dt);
}

error = EvolveOneTimestep(dt); // evolve the binary system one timestep

if (OPTIONS->EmitGravitationalRadiation()) {
EmitGravitationalWave(dt); // emit a graviataional wave
}

if (error != ERROR::NONE) { // SSE error for either constituent star?
evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution
}
Expand Down Expand Up @@ -2795,28 +2899,6 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning
}

if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution?

m_Star1->UpdatePreviousTimestepDuration();
m_Star2->UpdatePreviousTimestepDuration();

if (usingProvidedTimesteps) { // user-provided timesteps
// get new timestep
// - don't quantise
// - don't apply timestep multiplier
// (we assume user wants the timesteps in the file)
dt = timesteps[stepNum];
}
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 (dt < NUCLEAR_MINIMUM_TIMESTEP) {
dt = NUCLEAR_MINIMUM_TIMESTEP; // but not less than minimum
}
stepNum++; // increment stepNum
}
}

if (usingProvidedTimesteps && timesteps.size() > stepNum) { // all user-defined timesteps consumed?
Expand Down
16 changes: 12 additions & 4 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,11 +379,14 @@ class BaseBinaryStar {

double m_TotalEnergy;

double m_OrbitalAngularMomentumPrev;
double m_OrbitalAngularMomentum;
double m_OrbitalAngularMomentumPrev;
double m_OrbitalAngularMomentum;

double m_OrbitalEnergyPrev;
double m_OrbitalEnergy;
double m_DaDtGW; // Change in semi-major axis per time due to gravitational radiation
double m_DeDtGW; // Change in eccentricity per time due to gravitational radiation

double m_OrbitalEnergyPrev;
double m_OrbitalEnergy;

double m_ZetaLobe;
double m_ZetaStar;
Expand Down Expand Up @@ -413,6 +416,11 @@ 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()); }

void CalculateGravitationalRadiation();
void EmitGravitationalWave(const double p_Dt);

double ChooseTimestep(double p_Dt);

void CalculateEnergyAndAngularMomentum();

double CalculateDEccentricityTidalDt(const DBL_DBL_DBL_DBL p_ImKlm, const BinaryConstituentStar* p_Star);
Expand Down
8 changes: 7 additions & 1 deletion src/Options.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,8 @@ void Options::OptionValues::Initialise() {
m_EvolveMainSequenceMergers = false;
m_EvolvePulsars = false;
m_EvolveUnboundSystems = true;

m_EmitGravitationalRadiation = false;

m_NatalKickForPPISN = false;

m_DetailedOutput = false;
Expand Down Expand Up @@ -814,6 +815,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(false),
("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),
Expand Down
Loading

0 comments on commit 6a9fa62

Please sign in to comment.