Skip to content

Commit

Permalink
Merge pull request #1093 from TeamCOMPAS/PickerUpdate
Browse files Browse the repository at this point in the history 
Update to Picker+ (2024) fits for Hirai & Mandel 2-stage CE; exposure of functionality needed for tidal computations
  • Loading branch information
@jeffriley
jeffriley authored Apr 19, 2024
2 parents e2c12cd + b94f7aa commit 5aa2d21
Show file tree
Hide file tree
Showing 19 changed files with 291 additions and 220 deletions.
17 changes: 11 additions & 6 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1583,18 +1583,23 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() {

// Two-stage common envelope, Hirai & Mandel (2022)
else if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) {
double convectiveEnvelopeMass1 = m_Star1->CalculateConvectiveEnvelopeMass();
double convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1;
std::tie(convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1) = m_Star1->CalculateConvectiveEnvelopeMass();
double radiativeIntershellMass1 = m_MassEnv1 - convectiveEnvelopeMass1;
double endOfFirstStageMass1 = m_Mass1Final + radiativeIntershellMass1;
double convectiveEnvelopeMass2 = m_Star2->CalculateConvectiveEnvelopeMass();
double convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2;
std::tie(convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2) = m_Star2->CalculateConvectiveEnvelopeMass();
double radiativeIntershellMass2 = m_MassEnv2 - convectiveEnvelopeMass2;
double endOfFirstStageMass2 = m_Mass2Final + radiativeIntershellMass2;

// Stage 1: convective envelope removal on a dynamical timescale; assumes lambda = 1.5, motivated by bottom panel of Figure 3 of Hirai & Mandel 2022, including internal energy
double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius();
double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius();
// Stage 1: convective envelope removal on a dynamical timescale; assumes lambda = lambda_He
double lambda1 = m_Star1->CalculateConvectiveEnvelopeLambdaPicker(convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1);
double lambda2 = m_Star1->CalculateConvectiveEnvelopeLambdaPicker(convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2);

double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda1 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius();
double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda2 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius();
double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; // assumes immediate circularisation at periastron at start of CE
double k4 = (endOfFirstStageMass1 * endOfFirstStageMass2);
double k4 = endOfFirstStageMass1 * endOfFirstStageMass2;
double aFinalRsol = k4 / (k1 + k2 + k3);
m_SemiMajorAxis = aFinalRsol * RSOL_TO_AU;

Expand Down
35 changes: 29 additions & 6 deletions src/BaseStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2208,7 +2208,7 @@ double BaseStar::CalculateMassLossRateRSGVinkSabhahit2023() const {
if (utils::Compare(logL, logLkink) < 0) {
logMdot = -8.0 + 0.7 * logL - 0.7 * logM;
}
else if (utils::Compare(logL, logLkink) >= 0) {
else {
logMdot = -24.0 + 4.77 * logL - 3.99 * logM;
}

Expand Down Expand Up @@ -3805,19 +3805,42 @@ void BaseStar::CalculateBindingEnergies(const double p_CoreMass, const double p_
* Calculate convective envelope binding energy for the two-stage Hirai & Mandel (2022) common envelope formalism
*
*
* double CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda)
* double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda)
*
* @param [IN] p_CoreMass Core mass of the star (Msol)
* @param [IN] p_TotalMass Total mass of the star (Msol)
* @param [IN] p_ConvectiveEnvelopeMass Mass of the convective outer envelope (Msol)
* @param [IN] p_Radius Radius of the star (Rsol)
* @param [IN] p_Lambda Dimensionless parameter defining the binding energy
* @param [IN] p_Lambda Lambda parameter for the convective envelope
* @return Binding energy (erg)
*/
double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda) {
return CalculateBindingEnergy(p_CoreMass, p_ConvectiveEnvelopeMass, p_Radius, p_Lambda);
double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda) {
return CalculateBindingEnergy(p_TotalMass - p_ConvectiveEnvelopeMass, p_ConvectiveEnvelopeMass, p_Radius, p_lambda);
}


/*
* Approximates the lambda binding energy parameter of the outer convective envelope
*
* This is needed for the Hirai & Mandel (2022) two-stage CE formalism.
* Follows the fits of Picker, Hirai, Mandel (2024), arXiv:2402.13180 for lambda_He
*
*
* double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const
*
* @param [IN] p_convectiveEnvelopeMass Mass of the outer convective envelope shell
* @param [IN] p_maxConvectiveEnvelopeMass Maximum mass of the outer convective envelope shell at the onset of carbon burning
* @return Lambda binding energy parameter for the outer convective envelope
*/
double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const {

double m2 = 0.0023 * m_Log10Metallicity * m_Log10Metallicity + 0.0088 * m_Log10Metallicity + 0.013; // Eq. (12) and Table 1 of Picker, Hirai, Mandel (2024)
double b1 = m2 * m_Mass - 0.23; // Eq. (11) of Picker+ (2024)
double logLambda = p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass > 0.3 ?
0.42 * p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass + b1 : 0.3 * 0.42 + b1;

return exp(logLambda);
}

///////////////////////////////////////////////////////////////////////////////////////
// //
// MISCELLANEOUS FUNCTIONS / CONTROL FUNCTIONS //
Expand Down
14 changes: 9 additions & 5 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,13 @@ class BaseStar {

void CalculateBindingEnergies(const double p_CoreMass, const double p_EnvMass, const double p_Radius);

double CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda);
virtual double CalculateConvectiveCoreMass () const {return 0.0;}
virtual double CalculateConvectiveCoreRadius () const {return 0.0;}

virtual double CalculateConvectiveEnvelopeMass() const { return 0.0; }
double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda);

double CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const;
virtual DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple<double, double> (0.0, 0.0); }

virtual double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate);
virtual double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const { return 0.0; } // Default is 0.0
Expand Down Expand Up @@ -202,6 +206,8 @@ class BaseStar {

double CalculateRadialExpansionTimescale() const { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables

virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // default for stars with no convective envelope

void CalculateSNAnomalies(const double p_Eccentricity);

double CalculateSNKickMagnitude(const double p_RemnantMass, const double p_EjectaMass, const STELLAR_TYPE p_StellarType);
Expand Down Expand Up @@ -428,7 +434,7 @@ class BaseStar {

virtual double CalculateCOCoreMassAtPhaseEnd() const { return m_COCoreMass; } // Default is NO-OP
virtual double CalculateCOCoreMassOnPhase() const { return m_COCoreMass; } // Default is NO-OP

virtual double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // Default is NO-OP
static double CalculateCoreMassGivenLuminosity_Static(const double p_Luminosity, const DBL_VECTOR &p_GBParams);
virtual double CalculateCoreMassOnPhase() const { return m_CoreMass; } // Default is NO-OP
Expand Down Expand Up @@ -534,8 +540,6 @@ class BaseStar {
const double p_RadiusPrev,
const double p_DtPrev);

virtual double CalculateRadialExtentConvectiveEnvelope() const { return m_Radius; } // default for stars with no convective envelope

virtual double CalculateRadiusAtPhaseEnd() const { return m_Radius; } // Default is NO-OP
double CalculateRadiusAtZAMS(const double p_MZAMS) const;
virtual double CalculateRadiusOnPhase() const { return m_Radius; } // Default is NO-OP
Expand Down
20 changes: 13 additions & 7 deletions src/BinaryConstituentStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,9 @@ void BinaryConstituentStar::SetPostCEEValues() {
* m_CEDetails.CoreMass
* m_CEDetails.bindingEnergy
* m_CEDetails.lambda
*
* m_CEDetails.convectiveEnvelopeMass
* m_CEDetails.radiativeIntershellMass
* m_CEDetails.convectiveEnvelopeBindingEnergy
*
* void CalculateCommonEnvelopeValues()
*/
Expand All @@ -247,10 +249,6 @@ void BinaryConstituentStar::CalculateCommonEnvelopeValues() {
m_CEDetails.CoreMass = CoreMass();

m_CEDetails.lambda = 0.0; // default

m_CEDetails.convectiveEnvelopeMass = CalculateConvectiveEnvelopeMass();
m_CEDetails.radiativeIntershellMass = Mass() - CoreMass() - m_CEDetails.convectiveEnvelopeMass;
m_CEDetails.convectiveEnvelopeBindingEnergy = 0.0;

switch (OPTIONS->CommonEnvelopeLambdaPrescription()) { // which common envelope lambda prescription?

Expand Down Expand Up @@ -278,15 +276,23 @@ void BinaryConstituentStar::CalculateCommonEnvelopeValues() {
m_CEDetails.lambda = Lambda_Dewi();
m_CEDetails.bindingEnergy = BindingEnergy_Dewi();
break;

default: // unknown prescription
SHOW_WARN(ERROR::UNKNOWN_CE_LAMBDA_PRESCRIPTION, "Lambda = 0.0"); // show warning
}

if (utils::Compare(m_CEDetails.lambda, 0.0) <= 0) m_CEDetails.lambda = 0.0; // force non-positive lambda to 0

m_CEDetails.lambda *= OPTIONS->CommonEnvelopeLambdaMultiplier(); // multiply by constant (program option, default = 1.0)
m_CEDetails.convectiveEnvelopeBindingEnergy = CalculateConvectiveEnvelopeBindingEnergy(CoreMass(), m_CEDetails.convectiveEnvelopeMass, Radius(), m_CEDetails.lambda);

// Properties relevant for the Hirai & Mandel (2022) formalism
double maxConvectiveEnvelopeMass;
std::tie(m_CEDetails.convectiveEnvelopeMass, maxConvectiveEnvelopeMass) = CalculateConvectiveEnvelopeMass();
m_CEDetails.radiativeIntershellMass = Mass() - CoreMass() - m_CEDetails.convectiveEnvelopeMass;
if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE )
m_CEDetails.lambda = CalculateConvectiveEnvelopeLambdaPicker(m_CEDetails.convectiveEnvelopeMass, maxConvectiveEnvelopeMass);

m_CEDetails.convectiveEnvelopeBindingEnergy = CalculateConvectiveEnvelopeBindingEnergy(Mass(), m_CEDetails.convectiveEnvelopeMass, Radius(), m_CEDetails.lambda);
}


Expand Down
3 changes: 1 addition & 2 deletions src/CHeB.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ class CHeB: virtual public BaseStar, public FGB {

double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(CHeB) = 0.0

double CalculateConvectiveCoreRadius () const { return CalculateRemnantRadius (); } // Last paragraph of section 6 of Hurley+ 2000
double CalculateCoreMassAtPhaseEnd() const { return CalculateCoreMassAtBAGB(m_Mass0); } // Use class member variables
double CalculateCoreMassOnPhase(const double p_Mass, const double p_Tau) const;
double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Tau); } // Use class member variables
Expand All @@ -85,8 +86,6 @@ class CHeB: virtual public BaseStar, public FGB {
double CalculateLuminosityOnPhase(const double p_Mass, const double p_Tau) const;
double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Mass0, m_Tau); }

double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // CHeB stars don't have a convective envelope

double CalculateRadiusAtBluePhaseEnd(const double p_Mass) const;
double CalculateRadiusAtBluePhaseStart(const double p_Mass) const;

Expand Down
4 changes: 2 additions & 2 deletions src/FGB.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,9 +56,9 @@ class FGB: virtual public BaseStar, public HG {
double CalculateLuminosityOnPhase(const double p_Time) const;
double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age); } // Use class member variables

double CalculateRadialExtentConvectiveEnvelope() const { return GiantBranch::CalculateRadialExtentConvectiveEnvelope(); } // Skip HG
double CalculateRadialExtentConvectiveEnvelope() const { return m_Radius - CalculateConvectiveCoreRadius(); } // Skip HG

double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG - same as on phase
double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG - same as on phase
double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass, m_Luminosity); } // Use class member variables
double CalculateRadiusOnPhase(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG
double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass, m_Luminosity); } // Use class member variables
Expand Down
Loading

0 comments on commit 5aa2d21

Please sign in to comment.