Implement Coen's new mass transfer stability prescriptions

[SimonStevenson]

Is your feature request related to a problem? Please describe. We should implement Coen's new mass transfer stability prescription into COMPAS.

Describe the solution you'd like This consists of a few things: (@cneijssel correct me if I miss anything)

1. Allow user to change the definition of convective/radiative envelope from the current prescription (based on stellar type) to one based on effective temperature.
2. Allow user to change from using fixed zeta parameters for different stellar types to one based on results from Ge et al. 2015. This could either i) interpolate the value as a function of position on the HR diagram from a table of stellar models, or ii) use an analytic fit to such models for zeta as a function of position on the HR diagram.
3. Update documentation accordingly.

Describe alternatives you've considered NA

Additional context See proj_mt_stability channel on Slack for more details.

[SimonStevenson]

Regarding the first point (convective vs radiative envelopes).

This requires adding two new program options.

Firstly, one to choose the prescription, something like convective_envelope_prescription, with two options, something like STELLAR_TYPE and EFFECTIVE_TEMPERATURE.

Secondly, one to set the threshold effective temperature for the new prescription. We can set the default to whatever value @cneijssel is using.

For the new prescription, presumably this should only apply for evolved stars (stellar types 2--6 inclusive , right? Low mass stars can be much hotter and still have convective envelopes.

[SimonStevenson]

Determining the envelope type is currently done on a star by star basis (by which I mean in e.g. HG.h, ms_gt_07.h, etc), in the function DetermineEnvelopeType.

[SimonStevenson]

I have started working on this on a branch on my fork https://github.com/SimonStevenson/COMPAS/tree/new_mt_stability

So far I've added the program options for 1) and have started implementing the new convective_envelope_prescription for the different stars.

Reminder to self: check units in temperature comparison

[SimonStevenson]

Regarding point 2:

I'm currently a bit confused about our current zeta options, may need to ping @avigna and @cneijssel about this.

We currently have a common_envelope_zeta_prescription (which I believe is misnamed and should really be mass_transfer_zeta_prescription).

Somehow, even though our default choice of "common_envelope_zeta_prescription" is SOBERMAN, it seems that this only applies to Giants (so should be called something more related to giants).

In addition, it seems that the choices of zetaHertzsprungGap and zetaMainSequence are hardcoded, seperately from the choice of "common_envelope_zeta_prescription"

I was hoping that we could just add a new option "Ge2015" to "common_envelope_zeta_prescription", but I don't think this will work currently.

It looks like it might be better to rebuild the whole way mass transfer stability is determined in COMPAS.

Thoughts?

[SimonStevenson]

Most of the mess is in BaseBinaryStar::CalculateMassTransfer

[avigna]

Hi Simon, most of what you say is alright. Some comments (IMO): 1) Agreed about convective/radiatively envelope should just be a property of the star and not have to do with stability. Can be used for stability if needed. 2) Ge zetas were hardcoded. I remember doing minor tests for just defining zetas and q's and failing. Given my time constrains I ended up hardcoding them... but there should be a list of zetas, e.g. by stellar type, that one could define. Similar for mass ratio (I recall you started working on this at some point). 2.1) If we implement separate zetas, then Ge2015 would just pick a value for them. 3) Soberman is used for giants and naked helium stars (type = 3,4,5,6,7). It might be used for type =8,9 if force_case_BB_BC_stability=False, but not 100% sure. 4) At the end, there should be a function which returns true of false for stable mass transfer. The former leads to the mass transfer function, the latter to the CEE function.

(Probably more comments to come in the future.)

[cneijssel]

Zetas, stability, beta, gamma and the way we solve the orbit was done in multiple places with multiple projects/ideas behind it. Not familiar with the status of mass transfer and the stability in the new code. But given that it is ported from Legacy, it might have taken with it the entangled state of things. The naming of parameters made sense to me given the past prescriptions that we did. However, if we now want to consistently disentangle then it will be quite some work.

Not sure what you had in mind for the amount of work, but I guess from scratch would be best.

I have python code that goes through all the steps, although it does this for post-processing and could be cleaned up and used as a blueprint for possible restructure. But I understand we all have different ideas. This is one of the reasons why I started doing my project outside COMPAS-legacy

I think if we want to restructure / clean up we first create a pseudo code so we can look at dependencies of parameters. Writing that is straightforward since we all know the steps, allthough we might prefer different namings and order. Anyway my pseudo code is something along.

But implementing this will require hours of input in the code.

##########pseudo code order main parameters##

beta = determineBeta(stellarTypes, stellarTimeScales, programOptions) gamma = determineGamma(binaryOrbit, stellarMasses, programOptions) zeta_ad = determineZetaAdiabaticStar(beta, gamma, stellarTypes, stellarTemperatures, programOptions) zeta_RL = determineZetaRocheLobe(q, beta, gamma)

then determine stability and solve accordingly

pseudo code for determinZetaAdiabatic

determineZetaAdiabatic(beta, gamma, DonorParameters, programOptions){ ....* determine zeta adiabatic based on prescription*\

...if (programOptions.zetaAd = "stellarType_Fixed"){ ......if donorType MS { zeta_ad = programOptions.MSzeta} ......elif donorType HG{zeta_ad = programOptions.HGzeta} ......elif donorType FGB<TPAGB {zeta_ad = calculateZetaAdSoberMan(M, Mc)} ...}

...if (programOptions.zetaAd = "Temperature_Fixed"){ .........if donorTemperature >= logTconvective:{ ................if donorType MS { zeta_ad = programOptions.MSzeta} ................else donorType HG{zeta_ad = programOptions.HGzeta} .........} .........else{zeta_ad = calculateZetaAdSoberMan(M, Mc) .........} ...}

...if (programOptions.zetaAd = "Temperature_Ge"){ .........if donorTemperature >= logTconvective:{ ................zeta_ad = calculateZetaAdGe(R, MZAMS) .........} .........else{zeta_ad = calculateZetaAdSoberMan(M, Mc) .........} ...}

[cneijssel]

Sorry if I repeat some of Alejandro's remarks, forgot to press enter, had the text ready some time ago.

[avigna]

One of the things we should try to do to generalize (thinking forward) is getting away from stellar types (because of future addition of different stellar tracks) and quantities that rely on ZAMS (because of mass transfer). A least elegant approach is to do define stability parameters (zetas, q's) per stellar grid, but there might be a more elegant solution. I don't think HRD position in itself is enough.

[SimonStevenson]

@avigna Response to first post

1. Agreed and done.
2. Yeah I think there are a few choices we want. i) Arbitrary zetas for different stellar types ii) Zetas from models (e.g. SOBERMAN) iii) Zetas interpolated from data files (Ge2015) iv) Arbitrary critical mass ratios for different stellar types
3. Thanks for the clarification. My issue is with the naming of the option and the lack of flexibility; I will fix this in my restructuring (see also 2)
4. Yes, the idea would be to have a bool DetermineMassTransferStability() function which does exactly that (or possibly we could define a variable type STABILITY, which can take values STABLE and UNSTABLE, in the same way we have ENVELOPE which can be RADIATIVE or CONVECTIVE)

[SimonStevenson]

@cneijssel Thanks for the details. At the moment, mass transfer in COMPAS is the same as it always was. I want to generalise it and make it both more readable and more flexible whilst incorporating the option of using your new prescription. Your pseudo code should be very helpful for that, so thanks!

[SimonStevenson]

@avigna Response to second post: I think it will be difficult at the moment to have something which relies on anything other than basic stellar properties (mass, HRD diagram position) since that's all we really have. Even if we compute a grid of stellar models with zetas, they still won't account for binary interactions, so I don't think this is really feasible at the moment. However I don't think that this is a big issue, as I think HRD position should capture at least 1st order binary interaction effects (e.g. accreting mass moves you to a more massive star track, increases the luminosity of the star, changes effective zeta).

[SimonStevenson]

Here is some pseudo-code for the new structure of the mass transfer function:

We first check if we are allowing mass transfer, and then check if only one star is filling its Roche Lobe

Firstly, we need to determine if the mass transfer is stable

// Wrapper function to determine stability
STABILITY DetermineMassTransferStability(){

    // Create a variable to hold result, default is STABLE
    STABILITY stability = STABILITY::STABLE;

    switch (OPTIONS->MassTransferStabilityPrescription()) {

        case MASS_TRANSFER_STABILITY_PRESCRIPTION::NAME_OF_PRESCRIPTION:

            stability = DetermineMassTransferStabilityNAME_OF_PRESCRIPTION();

        // Add your new favourite mass transfer stability prescription here

        default : 
            // Show warning
            SHOW_WARN(ERROR::UNKNOWN_MASS_TRANSFER_STABILITY_PRESCRIPTION, "Using stability = STABLE");

    return stability;

}

// Then call this function as follows
STABILITY stability = DetermineMassTransferStability();

// Example of a specific function that determines mass transfer
// stability for a specific option (prescription)
STABILITY DetermineMassTransferStabilityCriticalMassRatios(){

    // Create a variable to hold result, default is STABLE
    STABILITY stability = STABILITY::STABLE;

    // Code that tests for stellar types and compares them to
    // user specified mass ratios goes here. Updates stability

    return stability;

}

Then once we have determined the stability, we decide what to do (similar to what Alejandro said)

switch (stability){
    case STABILITY::STABLE:
        DoStableMassTransfer();
    case STABILITY::UNSTABLE:
        GoToCommonEnvelopeEvolution();
}

For the moment, I will leave the rest of the way mass transfer is handelled as it is (the DoStableMassTransfer function above). In the future we should probably revisit that as well, but I don't want to bite off more than I can chew here.

So then the question is just what are the MassTransferStabilityPrescriptions that we want to implement (should include what we already have + Coen's new one).

Currently there are 4 options for the ill-named CommonEnvelopeZetaPrescription

{ CE_ZETA_PRESCRIPTION::STARTRACK, "STARTRACK" },
{ CE_ZETA_PRESCRIPTION::SOBERMAN,  "SOBERMAN" },
{ CE_ZETA_PRESCRIPTION::HURLEY,    "HURLEY" },
{ CE_ZETA_PRESCRIPTION::ARBITRARY, "ARBITRARY" }

These are implemented in GiantBranch::CalculateZeta

I'm not really sure what some of these do, but I can check. It seems that the STARTRACK method is not even implemented, so I will just remove that completely.

But then there are additional flags to force caseBB mass transfer to always be stable

--alwaysStableCaseBBBCFlag [=arg(=1)] (=0) Choose case BB/BC mass transfer to be always stable (default = FALSE) --forceCaseBBBCStabilityFlag [=arg(=1)] (=0) Force case BB/BC mass transfer to be only stable or unstable (default = FALSE)

and these fixed zetas

--zeta-adiabatic-arbitrary arg (=10000) Value of mass-radius exponent zeta adiabatic (default = 10000.000000) --zeta-thermal-arbitrary arg (=10000) Value of mass-radius exponent zeta adiabatic (default = 10000.000000) --zeta-hertzsprung-gap arg (=6.5) Value of mass-radius exponent zeta on the hertzstrpung gap (default = 6.500000) --zeta-main-sequence arg (=2) Value of mass-radius exponent zeta on the main sequence (default = 2.000000)

I think one important question is whether we want the prescription to be for all stellar types, or do we want to be able to select individually for each stellar type. E.g. do we want to be able to use SOBERMAN zetas for giants but a critical mass ratio for main sequence stars? That will require more options, but fundamentally can be achieved in a similar way.

Then we can implement

DetermineMassTransferStabilityNAME_OF_PRESCRIPTIONS

for each of these. How does this sound @avigna @cneijssel @jeffriley ?

[SimonStevenson]

Looks like what I refer to as STABILITY above already exists as MASS_TRANSFER

// Mass Transfer types
enum class MASS_TRANSFER: int { NONE, STABLE_TIMESCALE_NUCLEAR, STABLE_TIMESCALE_THERMAL, UNSTABLE_TIMESCALE_THERMAL, UNSTABLE_TIMESCALE_DYNAMICAL};
const COMPASUnorderedMap<MASS_TRANSFER, std::string> MASS_TRANSFER_LABEL = {
    { MASS_TRANSFER::NONE,                         "No mass transfer" },
    { MASS_TRANSFER::STABLE_TIMESCALE_NUCLEAR,     "Nuclear timescale stable mass transfer" },
    { MASS_TRANSFER::STABLE_TIMESCALE_THERMAL,     "Thermal timescale stable mass transfer" },
    { MASS_TRANSFER::UNSTABLE_TIMESCALE_THERMAL,   "Thermal timescale unstable mass transfer" },
    { MASS_TRANSFER::UNSTABLE_TIMESCALE_DYNAMICAL, "Dynamical timescale unstable mass transfer" }
};

[SimonStevenson]

I don't see where these are used though ...

[jeffriley]

They're not as far as I can tell. They were in constants.h in the legacy code (and only there I think), and I probably copied them over before I realised they weren't actually used. I probably left them there because I thought they might be useful at some stage (but more probably through a reluctance to accept I had wasted some time copying something that wasn't useful...).

[SimonStevenson]

Thanks Jeff! I think they will be useful! I'd appreciate your thoughts on the rest of this issue too when you have the time :)

[jeffriley]

Simon's pseudo-code looks like a good plan to me - we can refine if necessary as we see the actual implementation.

[ilyamandel]

Sorry, I didn't notice this issue earlier, so inadvertently did some similar work to @SimonStevenson 's fork. I added the FIXED_TEMPERATURE prescription to DetermineEnvelopeType() in version 02.12.05. I still need to test the results carefully (got distracted by other code changes since) and add the Ge+ interpolated zeta calculations. So keeping this issue open until this is checked.

[ilyamandel]

@SimonStevenson / @reinhold-willcox / @themikelau -- what's the status of this? I think much of this has been implemented already, right?

[reinhold-willcox]

@ilyamandel from a quick skim, I think we've implemented all that we need from the above discussion. Ge et al are nearly merged in from my other PR #849, and the zeta-based prescriptions generally have been superseded so I don't think we need to add more zeta prescriptions.

[ilyamandel]

OK, I'm pulling in #849 and closing this.