Skip to content

Commit a5794c1

Browse files
author
Ilya Mandel
committed
Fix to behaviour for stars that do not qualify for ECSN, see #1305 
1 parent bf16afa commit a5794c1

File tree

11 files changed

+81
-69
lines changed

11 files changed

+81
-69
lines changed

‎online-docs/COMPAS-2025methodsPaper.bib‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
@ARTICLE{COMPAS:2025,
2+
author = {{Team COMPAS: Mandel}, Ilya and {Riley}, Jeff and {Boesky}, Adam and {Brcek}, Adam and {Hirai}, Ryosuke and {Kapil}, Veome and {Lau}, Mike Y.~M. and {Merritt}, JD and {Rodr{\'\i}guez-Segovia}, Nicol{\'a}s and {Romero-Shaw}, Isobel and {Song}, Yuzhe and {Stevenson}, Simon and {Vajpeyi}, Avi and {van Son}, L.~A.~C. and {Vigna-G{\'o}mez}, Alejandro and {Willcox}, Reinhold},
3+
title = "{Rapid stellar and binary population synthesis with COMPAS: methods paper II}",
4+
journal = {arXiv e-prints},
5+
keywords = {Solar and Stellar Astrophysics, High Energy Astrophysical Phenomena, Instrumentation and Methods for Astrophysics},
6+
year = 2025,
7+
month = jun,
8+
eid = {arXiv:2506.02316},
9+
pages = {arXiv:2506.02316},
10+
doi = {10.48550/arXiv.2506.02316},
11+
archivePrefix = {arXiv},
12+
eprint = {2506.02316},
13+
primaryClass = {astro-ph.SR},
14+
adsurl = {https://ui.adsabs.harvard.edu/abs/2025arXiv250602316M},
15+
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
16+
}

‎online-docs/pages/how-to-cite.rst‎

Lines changed: 36 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,33 +1,36 @@
1-
Citing COMPAS
2-
-------------
3-
4-
If you use this code or parts of this code for results presented in a scientific publication, we would greatly appreciate you sending
5-
us your paper reference and making your input settings and output data publicly available by uploading it to the COMPAS Zenodo community.
6-
7-
Please also cite:
8-
9-
.. _cite-compas:
10-
11-
Team COMPAS: J. Riley `et al.` [:cite:year:`compas2021`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2021methodsPaper.bib>`
12-
13-
|br|
14-
We would also greatly appreciate an acknowledgement of the form:
15-
16-
"Simulations in this paper made use of the COMPAS rapid binary population synthesis code (version x.y.z), which is freely available at
17-
http://github.com/TeamCOMPAS/COMPAS."
18-
19-
|br|
20-
Furthermore,
21-
22-
If using the COMPAS model of gravitational wave selection effects, please cite :cite:t:`Barrett2018`.
23-
24-
If you use COMPAS's importance sampling algorithm STROOPWAFEL, please cite :cite:t:`Broekgaarden2019`.
25-
26-
If using COMPAS's integration over cosmic star formation history, please cite :cite:t:`Neijssel2019`.
27-
28-
If using the COMPAS model of (pulsational) pair instability supernova, please cite :cite:t:`Stevenson2019`.
29-
30-
If evolving pulsar spins and magnetic fields with COMPAS, please cite :cite:t:`Chattopadhyay2020`.
31-
32-
If you use the COMPAS model of chemically homogeneous evolution, please cite :cite:t:`Riley2021`.
33-
1+
Citing COMPAS
2+
-------------
3+
4+
If you use this code or parts of this code for results presented in a scientific publication, we would greatly appreciate you sending
5+
us your paper reference and making your input settings and output data publicly available by uploading it to the COMPAS Zenodo community.
6+
7+
Please also cite:
8+
9+
.. _cite-compas:
10+
11+
Team COMPAS: J. Riley `et al.` [:cite:year:`compas2021`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2021methodsPaper.bib>`
12+
13+
Team COMPAS: I. Mandel `et al.` [:cite:year:`compas2025`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2025methodsPaper.bib>`
14+
15+
16+
|br|
17+
We would also greatly appreciate an acknowledgement of the form:
18+
19+
"Simulations in this paper made use of the COMPAS rapid binary population synthesis code (version x.y.z), which is freely available at
20+
http://github.com/TeamCOMPAS/COMPAS."
21+
22+
|br|
23+
Furthermore,
24+
25+
If using the COMPAS model of gravitational wave selection effects, please cite :cite:t:`Barrett2018`.
26+
27+
If you use COMPAS's importance sampling algorithm STROOPWAFEL, please cite :cite:t:`Broekgaarden2019`.
28+
29+
If using COMPAS's integration over cosmic star formation history, please cite :cite:t:`Neijssel2019`.
30+
31+
If using the COMPAS model of (pulsational) pair instability supernova, please cite :cite:t:`Stevenson2019`.
32+
33+
If evolving pulsar spins and magnetic fields with COMPAS, please cite :cite:t:`Chattopadhyay2020`.
34+
35+
If you use the COMPAS model of chemically homogeneous evolution, please cite :cite:t:`Riley2021`.
36+

‎online-docs/pages/quick-links.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,4 +6,6 @@ Quick Links
66

77
./User guide/Program options/program-options-list-defaults
88
./User guide/COMPAS output/standard-logfiles-record-types
9+
./User guide/COMPAS output/standard-logfiles-record-specification-stellar.html
10+
./User guide/COMPAS output/standard-logfiles-record-specification-binary.html
911
./how-to-cite

‎src/EAGB.cpp‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1078,14 +1078,12 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_Force) {
10781078
*/
10791079
bool EAGB::ShouldSkipPhase() const {
10801080
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function
1081-
#define gbParams(x) m_GBParams[static_cast<int>(GBP::x)] // for convenience and readability - undefined at end of function
10821081

10831082
double McCOBAGB = CalculateCOCoreMassOnPhase(timescales(tHeI) + timescales(tHe));
1084-
double McSN = std::max(gbParams(McSN), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya**
1083+
double McSN = std::max(CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast<int>(GBP::McBAGB)]), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper
10851084

10861085
return (utils::Compare(McSN, m_COCoreMass) < 0); // skip phase if core is heavy enough to go supernova
10871086

1088-
#undef gbParams
10891087
#undef timescales
10901088
}
10911089

@@ -1121,14 +1119,12 @@ bool EAGB::ShouldEvolveOnPhase() const {
11211119
*/
11221120
bool EAGB::IsSupernova() const {
11231121
#define timescales(x) m_Timescales[static_cast<int>(TIMESCALE::x)] // for convenience and readability - undefined at end of function
1124-
#define gbParams(x) m_GBParams[static_cast<int>(GBP::x)] // for convenience and readability - undefined at end of function
11251122

11261123
double McCOBAGB = CalculateCOCoreMassOnPhase(timescales(tHeI) + timescales(tHe));
1127-
double McSN = std::max(gbParams(McSN), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya**
1124+
double McSN = std::max(CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast<int>(GBP::McBAGB)]), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya**
11281125

11291126
return (utils::Compare(McSN, m_COCoreMass) <= 0); // core is heavy enough to go Supernova
11301127

1131-
#undef gbParams
11321128
#undef timescales
11331129
}
11341130

‎src/GiantBranch.cpp‎

Lines changed: 5 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -276,9 +276,6 @@ void GiantBranch::CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams)
276276
gbParams(McBAGB) = CalculateCoreMassAtBAGB(p_Mass);
277277
gbParams(McDU) = CalculateCoreMassAt2ndDredgeUp_Static(gbParams(McBAGB));
278278
gbParams(McBGB) = CalculateCoreMassAtBGB(p_Mass, p_GBParams);
279-
280-
gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB));
281-
282279
#undef gbParams
283280
}
284281

@@ -341,8 +338,6 @@ void GiantBranch::CalculateGBParams_Static(const double p_Mass,
341338
gbParams(McBAGB) = CalculateCoreMassAtBAGB_Static(p_Mass, p_BnCoefficients);
342339
gbParams(McBGB) = CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams);
343340

344-
gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB));
345-
346341
#undef gbParams
347342
}
348343

@@ -813,16 +808,17 @@ double GiantBranch::CalculateCoreMassAtBGB_Static(const double p_Mass,
813808
/*
814809
* Calculate the core mass at which the Asymptotic Giant Branch phase is terminated in a SN/loss of envelope
815810
*
816-
* Hurley et al. 2000, eq 75 -- but note we use MECS rather than MCH
811+
* Hurley et al. 2000, eq 75 -- but note we may use MECS rather than MCH as the CO core mass threshold
817812
*
818813
*
819-
* double CalculateCoreMassAtSupernova_Static(const double p_McBAGB)
814+
* double CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB)
820815
*
816+
* @param [IN] p_Mthreshold Threshold mass, typically either MECS for stars that may explode in ECSNe or MCH otherwise
821817
* @param [IN] p_McBAGB Core mass at the Base of the Asymptotic Giant Branch in Msol
822818
* @return Maximum core mass before supernova on the Asymptotic Giant Branch
823819
*/
824-
double GiantBranch::CalculateCoreMassAtSupernova_Static(const double p_McBAGB) {
825-
return std::max(MECS, (0.773 * p_McBAGB) - 0.35);
820+
double GiantBranch::CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB) {
821+
return std::max(p_Mthreshold, (0.773 * p_McBAGB) - 0.35);
826822
}
827823

828824

‎src/GiantBranch.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
3636
double CalculateCoreMassAtBGB(const double p_Mass, const DBL_VECTOR &p_GBParams);
3737
static double CalculateCoreMassAtBGB_Static(const double p_Mass, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_GBParams);
3838
double CalculateCoreMassAtHeIgnition(const double p_Mass) const;
39-
static double CalculateCoreMassAtSupernova_Static(const double p_McBAGB);
39+
static double CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB);
4040

4141
static double CalculateCoreMass_Luminosity_B_Static(const double p_Mass);
4242
static double CalculateCoreMass_Luminosity_D_Static(const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs);

‎src/HeHG.cpp‎

Lines changed: 3 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -65,8 +65,6 @@ void HeHG::CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) {
6565
gbParams(McBAGB) = CalculateCoreMassAtBAGB();
6666
gbParams(McBGB) = CalculateCoreMassAtBGB(p_Mass, p_GBParams);
6767

68-
gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB));
69-
7068
#undef gbParams
7169
}
7270

@@ -134,8 +132,6 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0,
134132
gbParams(McBAGB) = p_Mass0;
135133
gbParams(McBGB) = GiantBranch::CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams);
136134

137-
gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB));
138-
139135
#undef gbParams
140136
}
141137

@@ -412,12 +408,10 @@ double HeHG::ChooseTimestep(const double p_Time) const {
412408
* @return Boolean flag: true if evolution on this phase should continue, false if not
413409
*/
414410
bool HeHG::ShouldEvolveOnPhase() const {
415-
#define gbParams(x) m_GBParams[static_cast<int>(GBP::x)] // for convenience and readability - undefined at end of function
416411

417412
double McMax = CalculateMaximumCoreMass(m_Mass);
418-
return ((utils::Compare(m_COCoreMass, McMax) <= 0 || utils::Compare(McMax, gbParams(McSN)) >= 0) && !ShouldEnvelopeBeExpelledByPulsations()); // Evolve on HeHG phase if McCO <= McMax or McMax >= McSN and envelope is not ejected by pulsations
419-
420-
#undef gbParams
413+
double McSN = CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast<int>(GBP::McBAGB)]);
414+
return ((utils::Compare(m_COCoreMass, McMax) <= 0 || utils::Compare(McMax, McSN) >= 0) && !ShouldEnvelopeBeExpelledByPulsations()); // Evolve on HeHG phase if McCO <= McMax or McMax >= McSN and envelope is not ejected by pulsations
421415
}
422416

423417

@@ -484,7 +478,7 @@ bool HeHG::IsSupernova() const {
484478
return (utils::Compare(m_Mass, MECS) > 0);
485479
}
486480

487-
return (utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(m_GBParams[static_cast<int>(GBP::McBAGB)])) >= 0); // Go supernova if CO core mass large enough
481+
return (utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast<int>(GBP::McBAGB)])) >= 0); // Go supernova if CO core mass large enough
488482
}
489483

490484
/*

‎src/TPAGB.cpp‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -978,8 +978,11 @@ double TPAGB::ChooseTimestep(const double p_Time) const {
978978
* @return Boolean flag: true if star has gone Supernova, false if not
979979
*/
980980
bool TPAGB::IsSupernova() const {
981-
// no supernova if CO core mass is too low or helium core mass is too low at base of AGB or the envelope has already been removed
982-
return utils::Compare(m_COCoreMass, m_GBParams[static_cast<int>(GBP::McSN)]) >= 0 &&
983-
utils::Compare(CalculateInitialSupernovaMass(), OPTIONS->MCBUR1()) >= 0 &&
984-
utils::Compare(m_COCoreMass, m_Mass) < 0;
981+
double snMass = CalculateInitialSupernovaMass();
982+
bool isCCSN = utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast<int>(GBP::McBAGB)])) >= 0 &&
983+
utils::Compare(snMass, OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0;
984+
bool isECSN = utils::Compare(snMass, MCBUR2) < 0 && (!m_MassTransferDonorHistory.empty() || OPTIONS->AllowNonStrippedECSN()) &&
985+
utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast<int>(GBP::McBAGB)])) >= 0 &&
986+
utils::Compare(snMass, OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0;
987+
return isCCSN || isECSN;
985988
}

‎src/TPAGB.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ class TPAGB: virtual public BaseStar, public EAGB {
4848

4949
// member functions - alphabetically
5050
DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple<double, double> (m_Mass-m_CoreMass, m_Mass-m_CoreMass); } // assume entire envelope is convective for TPAGB stars
51-
double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast<int>(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; }
51+
double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; }
5252
double CalculateCOCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // McCO(TPAGB) = Mc(TPAGB)Same as on phase
5353

5454
double CalculateConvectiveCoreRadius() const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000
@@ -103,7 +103,7 @@ class TPAGB: virtual public BaseStar, public EAGB {
103103
void ResolveHeliumFlash() { } // NO-OP
104104
STELLAR_TYPE ResolveSkippedPhase() { return m_StellarType; } // NO-OP
105105

106-
bool ShouldEvolveOnPhase() const { return ((utils::Compare(m_COCoreMass, std::min(m_GBParams[static_cast<int>(GBP::McSN)], m_Mass)) < 0) && !ShouldEnvelopeBeExpelledByPulsations()); } // Evolve on TPAGB phase if envelope is not lost and not going supernova
106+
bool ShouldEvolveOnPhase() const { return (utils::Compare(m_COCoreMass, m_Mass) < 0 && !IsSupernova() && !ShouldEnvelopeBeExpelledByPulsations()); } // Evolve on TPAGB phase if envelope is not lost and not going supernova
107107
bool ShouldSkipPhase() const { return false; } // Never skip TPAGB phase
108108

109109
};

‎src/changelog.h‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -944,10 +944,10 @@
944944
// - Cleaned up stability check functions in BaseBinaryStar.cpp for clarity, and to allow for critical mass ratios to be checked correctly
945945
// 02.33.01 RTW - Sep 26, 2022 - Defect repair:
946946
// - Fixed interpolation of MACLEOD_LINEAR gamma for specific angular momentum. Previously interpolated on the gamma value, now interpolates in orbital separation
947-
// 02.33.02 IM - Nov 27, 2022 - Defect repair:
947+
// 02.33.02 IM - Nov 27, 2022 - Defect repair:
948948
// - Fixed ignored value of input radius when computing the thermal timescale, relevant if using Roche lobe radius instead (issue #853)
949949
// - Cleaned code and comments around the use of MT_THERMALLY_LIMITED_VARIATION::RADIUS_TO_ROCHELOBE vs. C_FACTOR (issue #850)
950-
// 02.34.00 IM - Nov 28, 2022 - Enhancement:
950+
// 02.34.00 IM - Nov 28, 2022 - Enhancement:
951951
// - Adding framework for Hirai & Mandel 2-stage common envelope formalism
952952
// (placeholders for now -- will have identical results to default version)
953953
// - Placed Dewi CE prescription on parity with others
@@ -1579,9 +1579,12 @@
15791579
// 03.20.02 IM - May 30, 2026 - Defect repair, enhancement:
15801580
// - Included unit conversion in WhiteDwarfs::CalculateEtaPTY()
15811581
// - All critical mass ratios now return the HURLEY_HJELLMING_WEBBINK_QCRIT_WD for white dwarfs and 0 for other remnant donors (only stable mass transfer) as fix for issue #1385
1582+
// 03.20.03 IM - June 18, 2026 - Defect repair, enhancement:
1583+
// - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02)
1584+
// - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static()
1585+
// - Removed McSN from GBParams, instead computed on the fly when needed
15821586

1583-
1584-
const std::string VERSION_STRING = "03.20.02";
1587+
const std::string VERSION_STRING = "03.20.03";
15851588

15861589

15871590
# endif // __changelog_h__

0 commit comments

Comments
 (0)