Skip to content

Commit

Permalink
Allow for evolution of pulsar formed from main sequence stellar merger
Browse files Browse the repository at this point in the history
  • Loading branch information
@SimonStevenson
SimonStevenson committed Jan 20, 2025
1 parent 36b45ad commit 378ba5b
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 40 deletions.
89 changes: 51 additions & 38 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2850,45 +2850,52 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {
*/
double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {

double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()); // timestep based on orbital timescale
// Timesteps required by individual stars
double dt1 = m_Star1->CalculateTimestep();
double dt2 = m_Star2->CalculateTimestep();
double dt = std::min(dt1, dt2);

if (OPTIONS->EmitGravitationalRadiation()) { // emitting GWs?
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 (!IsUnbound()){ // Check that binary is bound

if (OPTIONS->EmitGravitationalRadiation()) { // emitting GWs?
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::KAPIL2024) { // tides prescription = KAPIL2024
// yes - need to adjust dt

double omega = OrbitalAngularVelocity();

DBL_DBL_DBL_DBL ImKlm1 = m_Star1->CalculateImKlmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
DBL_DBL_DBL_DBL ImKlm2 = m_Star2->CalculateImKlmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());

double DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKlm1, m_Star1);
double DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKlm2, m_Star2);

double DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(ImKlm1, m_Star1);
double DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(ImKlm2, m_Star2);

double DOmega1Dt_tidal = CalculateDOmegaTidalDt(ImKlm1, m_Star1);
double DOmega2Dt_tidal = CalculateDOmegaTidalDt(ImKlm2, m_Star2);

// Ensure that the change in orbital and spin properties due to tides in a single timestep is constrained (to 1 percent by default)
// Limit the spin evolution of each star based on the orbital frequency rather than its spin frequency, since tides should not cause major problems until synchronization.
double Dt_SemiMajorAxis1_tidal = utils::Compare(DSemiMajorAxis1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis1Dt_tidal) * YEAR_TO_MYR;
double Dt_SemiMajorAxis2_tidal = utils::Compare(DSemiMajorAxis2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis2Dt_tidal) * YEAR_TO_MYR;
double Dt_SemiMajorAxis_tidal = std::min(Dt_SemiMajorAxis1_tidal, Dt_SemiMajorAxis2_tidal);

double Dt_Eccentricity1_tidal = utils::Compare(DEccentricity1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity1Dt_tidal) * YEAR_TO_MYR;
double Dt_Eccentricity2_tidal = utils::Compare(DEccentricity2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity2Dt_tidal) * YEAR_TO_MYR;
double Dt_Eccentricity_tidal = std::min(Dt_Eccentricity1_tidal, Dt_Eccentricity2_tidal);

double Dt_Omega1_tidal = utils::Compare(DOmega1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega1Dt_tidal) * YEAR_TO_MYR;
double Dt_Omega2_tidal = utils::Compare(DOmega2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega2Dt_tidal) * YEAR_TO_MYR;
double Dt_Omega_tidal = std::min(Dt_Omega1_tidal, Dt_Omega2_tidal);

dt = std::min(dt, std::min(Dt_SemiMajorAxis_tidal, std::min(Dt_Eccentricity_tidal, Dt_Omega_tidal)));
if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::KAPIL2024) { // tides prescription = KAPIL2024
// yes - need to adjust dt

double omega = OrbitalAngularVelocity();

DBL_DBL_DBL_DBL ImKlm1 = m_Star1->CalculateImKlmTidal(omega, m_SemiMajorAxis, m_Star2->Mass());
DBL_DBL_DBL_DBL ImKlm2 = m_Star2->CalculateImKlmTidal(omega, m_SemiMajorAxis, m_Star1->Mass());

double DSemiMajorAxis1Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKlm1, m_Star1);
double DSemiMajorAxis2Dt_tidal = CalculateDSemiMajorAxisTidalDt(ImKlm2, m_Star2);

double DEccentricity1Dt_tidal = CalculateDEccentricityTidalDt(ImKlm1, m_Star1);
double DEccentricity2Dt_tidal = CalculateDEccentricityTidalDt(ImKlm2, m_Star2);

double DOmega1Dt_tidal = CalculateDOmegaTidalDt(ImKlm1, m_Star1);
double DOmega2Dt_tidal = CalculateDOmegaTidalDt(ImKlm2, m_Star2);

// Ensure that the change in orbital and spin properties due to tides in a single timestep is constrained (to 1 percent by default)
// Limit the spin evolution of each star based on the orbital frequency rather than its spin frequency, since tides should not cause major problems until synchronization.
double Dt_SemiMajorAxis1_tidal = utils::Compare(DSemiMajorAxis1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis1Dt_tidal) * YEAR_TO_MYR;
double Dt_SemiMajorAxis2_tidal = utils::Compare(DSemiMajorAxis2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_SemiMajorAxis / DSemiMajorAxis2Dt_tidal) * YEAR_TO_MYR;
double Dt_SemiMajorAxis_tidal = std::min(Dt_SemiMajorAxis1_tidal, Dt_SemiMajorAxis2_tidal);

double Dt_Eccentricity1_tidal = utils::Compare(DEccentricity1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity1Dt_tidal) * YEAR_TO_MYR;
double Dt_Eccentricity2_tidal = utils::Compare(DEccentricity2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * m_Eccentricity / DEccentricity2Dt_tidal) * YEAR_TO_MYR;
double Dt_Eccentricity_tidal = std::min(Dt_Eccentricity1_tidal, Dt_Eccentricity2_tidal);

double Dt_Omega1_tidal = utils::Compare(DOmega1Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega1Dt_tidal) * YEAR_TO_MYR;
double Dt_Omega2_tidal = utils::Compare(DOmega2Dt_tidal, 0.0) == 0 ? dt : std::abs(TIDES_MAXIMUM_ORBITAL_CHANGE_FRAC * omega / DOmega2Dt_tidal) * YEAR_TO_MYR;
double Dt_Omega_tidal = std::min(Dt_Omega1_tidal, Dt_Omega2_tidal);

dt = std::min(dt, std::min(Dt_SemiMajorAxis_tidal, std::min(Dt_Eccentricity_tidal, Dt_Omega_tidal)));
}
}

dt *= p_Multiplier;

return std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, TIDES_MNIMUM_FRACTIONAL_NUCLEAR_TIME * NUCLEAR_MINIMUM_TIMESTEP); // quantised and not less than minimum
Expand Down Expand Up @@ -3199,8 +3206,14 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
else if (!OPTIONS->EvolveDoubleWhiteDwarfs() && IsWDandWD()) { // double WD and their evolution is not enabled?
evolutionStatus = EVOLUTION_STATUS::WD_WD; // yes - do not evolve double WD systems
}
else if ((HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !OPTIONS->EvolveMainSequenceMergers()) || IsMRandRemant()) { // at least one massless remnant and not evolving MS merger products, or is MR + stellar remnant
evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT; // yes - stop evolution
else if ((HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT }) && !OPTIONS->EvolveMainSequenceMergers()) ||
IsMRandRemant()) { // at least one massless remnant and not evolving MS merger products, or is MR + stellar remnant
if (IsMRandNS() && OPTIONS->EvolvePulsars()){ // However, keep evolving if the stellar remnant is a neutron star and we are evolving pulsars
evolutionStatus = EVOLUTION_STATUS::CONTINUE;
}
else{
evolutionStatus = EVOLUTION_STATUS::MASSLESS_REMNANT; // yes - stop evolution
}
}
}
}
Expand Down
1 change: 1 addition & 0 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,7 @@ class BaseBinaryStar {
}
bool IsNSandBH() const { return HasOneOf({STELLAR_TYPE::NEUTRON_STAR}) && HasOneOf({STELLAR_TYPE::BLACK_HOLE}); }
bool IsNSandNS() const { return HasTwoOf({STELLAR_TYPE::NEUTRON_STAR}); }
bool IsMRandNS() const { return HasOneOf({STELLAR_TYPE::NEUTRON_STAR}) && HasOneOf({STELLAR_TYPE::MASSLESS_REMNANT}); }
bool IsMRandRemant() const { return HasOneOf({STELLAR_TYPE::MASSLESS_REMNANT}) && HasOneOf({STELLAR_TYPE::HELIUM_WHITE_DWARF, STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF, STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF, STELLAR_TYPE::NEUTRON_STAR, STELLAR_TYPE::BLACK_HOLE}); }
bool IsUnbound() const { return (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || (utils::Compare(m_Eccentricity, 1.0) > 0)); } // semi major axis <= 0.0 means unbound, presumably by SN)
bool IsWDandWD() const { return HasTwoOf({STELLAR_TYPE::HELIUM_WHITE_DWARF, STELLAR_TYPE::CARBON_OXYGEN_WHITE_DWARF, STELLAR_TYPE::OXYGEN_NEON_WHITE_DWARF}); }
Expand Down
2 changes: 1 addition & 1 deletion src/MR.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ class MR: virtual public BaseStar, public Remnants {
double CalculateMomentOfInertia() const { return 0.0; } // No moment of inertia for massless remnants - use 0.0
double CalculateMomentOfInertiaAU() const { return 0.0; } // No moment of inertia for massless remnants - use 0.0

double CalculateTimestep() const { return std::numeric_limits<double>::max(); } // Can take arbitrarily long time steps for massless remnants -- nothing is happening
double ChooseTimestep(const double p_Time) const { return std::numeric_limits<double>::max(); } // Can take arbitrarily long time steps for massless remnants -- nothing is happening

void SetPulsarParameters() const { } // NO-OP

Expand Down
5 changes: 4 additions & 1 deletion src/changelog.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1443,7 +1443,10 @@
// - Added treatment for rejuvenation of main sequence accretors when the new prescription is used
// 03.12.01 JR - Jan 17, 2025 - Defect repair:
// - (partial?) fix for issue #1149 - remove conditional from TPAGB::IsSupernova(). Whether it fixes issue 1149 completely or not, the conditional shouldn't be there...
// 03.12.02 SS - Jan 20, 2025 - Defect repair:
// - fix for issue #1324 - allow for evolution of pulsars formed from main sequence merger products. Changed CalculateTimestep to ChooseTimestep in MR.h and added check for bound binary to BaseBinaryStar::ChooseTimestep

const std::string VERSION_STRING = "03.12.01";

const std::string VERSION_STRING = "03.12.02";

# endif // __changelog_h__

0 comments on commit 378ba5b

Please sign in to comment.