Two Stage CE: unusual metallicity dependent pileup in yields

[jmerritt1]

When using the Two Stage Common Envelope formalism, and plotting BBH yields as a function of metallicity, we are finding two strange pileups:

Naively, there shouldn't be any Z dependence in the two stage formalism, correct? Our expectation would be a general supression of DCOs across the plot as compared to the Energy formalism, but the sharp features seem suspicious.

urgency_low - This issue is not urgent

severity_minor - This is a minor bug with minimal impact

To Reproduce Steps to reproduce the behavior:

1. run with the following .yaml file, varying --common-envelope-formalism compasExpYaml.txt (as a .yaml file, not .txt)
2. note that m1-min is 10.0

3. seed is 0

   • COMPAS [e.g. v03.01.07]

[ilyamandel]

Very interesting! There are two places where metallicity comes in explicitly in common envelope modelling, other than via stellar radii (and therefore possibly radial responses to mass loss, i.e., zetas), winds, etc.:

• In the classical alpha-lambda formalism, Nanjing lambda fits explicitly depend on metallicity and have a sharp transition, though @themikelau 's fits smoothed that;
• In the Picker fits, there is a dependence of the convective envelope mass on metallicity.

Here, it looks like the latter is the more likely culprit; I'll aim to make some plots.

[jmerritt1]

Trying to find some features of the systems in the pileups, it's clear that the lower Z population are forming from CE events with larger initial separations than usual. (all DCOs from a run are plotted, those that underwent CE are in blue).

[jmerritt1]

and double core CEs in yellow

[ilyamandel]

Very interesting. There also seems to be a band at metallicites between 0.004 and 0.007.

Is a@ZAMS in units of AU?

Good pointer to checking specifically on double-core CE -- I haven't investigated that specifically in the 2-stage CE context.

[jmerritt1]

Is a@ZAMS in units of AU?

Yes, it's in AU.

[jeffriley]

As a check, I ran two runs of COMPAS v03.01.07 as follows:

./compas --random-seed 0 -n 100000 --metallicity-distribution loguniform --metallicity-min 0.0001 --metallicity-max 0.001
./compas --random-seed 0 -n 100000 --metallicity-distribution loguniform --metallicity-min 0.0001 --metallicity-max 0.001 --common-envelope-formalism two_stage

I binned BHBH, merging in a Hubble time, into 10 linear metallicity bins, each of width 0.0001. The results were:

Note this was only for one end (the left) of the original plot @jmerritt1 posted, but I'm not seeing a big increase in yield for the two-stage common envelope. Maybe highlighted/exacerbated by binning choice?

[jmerritt1]

We discussed a few key differences: I do have logbinning in Z for the yields, and some of the other settings are non-default (notably CHE:NONE, etc. in the attached .yaml file)

[ilyamandel]

I think there are probably several issues to think about here, including connections to issues #1171 (metallicity-sensitive behaviour) and #1180 (Picker fits). However, there is one thing that's definitely an issue, so maybe let's start there first:

How to handle double-core CEs, including in the 2-stage formalism? (Pinging @ryosuke-hirai , @SimonStevenson , @avigna )

As a quick reminder, we consider a mass transfer episode to be a double-core CE when at least one star overflows its Roche Lobe in a way that is determined to lead to dynamically unstable mass transfer and both stars are evolved (Hertzsprung Gap, giants, or helium HG or GB).

If that's the case, we remove both envelopes.

In the classical alpha-lambda formalism, the order of removal doesn't matter, and in fact we remove both envelopes at once: compute the sum of the binding energies and then change the orbit so the difference between the initial orbital energy and the orbital energy of the remaining cores matches the alpha-lambda prescription.

In the 2-stage prescription, we do exactly the same thing for the first stage (modulo replacing "envelopes" with "convective envelopes" and "cores" with "cores + radiative intershells" in the sentence above). But then, when dealing with thermal timescale mass transfer from the exposed radiative intershells, the order does matter. And at the moment, we arbitrarily transfer mass from star1 (the ZAMS primary) to star2 (the ZAMS secondary) first, then in the reverse. [There's a question of whether the mass of the primary is properly updated between these two mass transfers, but let's leave the implementation aside for the moment -- for now, I'd like to understand what the "right" thing to do is.]

I see two options:

1. In the second stage of a 2-stage double-core CE, transfer the radiative intershell mass first from the star that initially overflowed the Roche lobe, starting the entire CE process in the first place; then reduce that star's mass to its core mass before transferring the radiative intershell in the opposite direction, from the initial accretor onto this core. (Aside: if we go this route, I would likely take this opportunity to finally implement #273: Replace a double-core CE with two sequential single-core CEs, starting with the one from the original donor. We would then check if at the end of the CE, the secondary overflows its Roche Lobe. If it does and has an envelope, we immediate go into a second CE phase, this time originating from the secondary. If the secondary has no envelope, and either the secondary or the primary overflows the Roche lobe immediately after CE, we put up our usual RLOF-immediately-after-CE flag and treat this as a merger. Otherwise, the system survives CE without doing anything to the secondary’s tightly bound envelope.)

2. After stripping the outer convective envelopes, we attempt to estimate the thermal timescales of the remaining radiative intershells. The first one to be transferred is the one with the shorter thermal timescale. (Aside: this solution is not, I think, consistent with replacing a double-core CE with two sequential single-core CEs, so we would abandon #273 if we went this route.)

Thoughts?

[jeffriley]

Two more runs, this time with v03.01.09 (shouldn't make a difference here), this time with @jmerritt1's option values:

./compas --random-seed 0 -n 100000 --metallicity-distribution loguniform --metallicity-min 0.0001 --metallicity-max 0.001 --initial-mass-min 10.0 --allow-rlof-at-birth false --pulsational-pair-instability-prescription FARMER --chemically-homogeneous-evolution-mode None

./compas --random-seed 0 -n 100000 --metallicity-distribution loguniform --metallicity-min 0.0001 --metallicity-max 0.001 --initial-mass-min 10.0 --allow-rlof-at-birth false --pulsational-pair-instability-prescription FARMER --chemically-homogeneous-evolution-mode None --common-envelope-formalism two_stage

Same linear binning as above. Results:

So it looks like the jump in yields for the lower metallicities is there. Must e some interaction between the option values and the 2-stage CE formalism.

[ryosuke-hirai]

I don't have any strong opinions for double-core CEs, but option 1 sounds more straightforward to me. Trying to estimate the thermal timescales of each "core" introduces unnecessary complications where we don't even know what is supposed to happen.

As for the pileups, I also don't have any good explanation but it may be useful to look at T_onset/T_min as a function of metallicity and mass.

[avigna]

Hmn, tricky. In principle the double core should mostly work for binaries with mass ratio very close to unity. The picture I have is that a donor would rapidly transfer most of the convective envelope to the accretor, that would then inflate as a hamstar and lead to reverse RLOF. I assume this would happen rather quickly, on a similar timescale than the expansion of the radiative layer of the donor. This sounds complex, but maybe just initially removing the two envelopes as we do now makes sense.

If you depart from mass ratio of unity, then you might have like the two separate common envelopes. We could assume or estimate timescales... but these would be rather ad-hoc. We could also add some mass ratio in which we expect double-core CEE to occur... but this would also be ad-hoc.

So, no definitive answers, but hope these thoughts help.

[ilyamandel]

Thank you, @ryosuke-hirai and @avigna ! I will go with option 1 (I agree, it's simpler), hopefully in the next few days. And yes, investigating the Picker fits is also high on the coding priority list.
Still have 3 papers to finish reading though, the winds one being first, so thank you for your patience to all those waiting for me on code development!

[ilyamandel]

The option above has been implemented in #1226 .
@jmerritt1 -- this may not have resolved the issue you pointed out; feel free to check now or wait until I implement one more planned change, to the Picker fits.

[jmerritt1]

FYI, @jeffriley found that the pileups may be due to the effects of the Farmer PPI prescription in conjunction with 2-stage CE. Doing some comparison runs to confirm this.

[jeffriley]

@jeffriley found that the pileups may be due to the effects of the Farmer PPI prescription in conjunction with 2-stage CE.

I'm no longer convinced... I did some runs with COMPAS v03.01.09 before I left Aspen, but didn't get a chance to consolidate them until now. Here are my results (these are just raw yields of BBH that will merge in a Hubble time, for each of the 9 metallicity bins shown):

@jmerritt1 I don't think I see with my runs what you posted. Can you check the option values for your results, just to make sure they were consistent across your runs (if you didn't use a grid file, the option values listed in Run_Details are the ones that were used)?

EDIT: well, maybe there is some variation there - but I don't think I see the ~40% variation that @jmerritt1 posted. Bear in mind that these are linear metallicity bins (each 0.0001), whereas @jmerritt1 posted log bins.

[ilyamandel]

Thank you for checking this, @jeffriley !

Note that the changes in #1226 are relevant regardless: they may not fix this particular issue (which may not even be an issue, according to Jeff's results?), but they should improve how 2-stage CE envelopes are handled, so should probably become the starting point for future investigations.

[jeffriley]

Ok, so this is the plot Simon produced:

Using these commandlines:

./COMPAS -n 100000 --metallicity-distribution LOGUNIFORM --common-envelope-formalism ENERGY --initial-mass-min 10.0 --output-container DEFAULT > default.out &
./COMPAS -n 100000 --metallicity-distribution LOGUNIFORM --common-envelope-formalism TWO_STAGE --initial-mass-min 10.0 --output-container CHECK2CE > check2CEchanges.out &

I ran the same two runs with COMPAS v03.02.01 and got the following two plots:

I used the metallicities of all entries in the DCO file that had ST1 = BH, ST2 = BH, and Merges in a Hubble time TRUE.

The first uses 20 log bins over the range of metallicities; the second is 100 equal-size linear bins over the same range.

Since neither of us specified a starting random seed, our runs will vary - but that shouldn't be by too much.

I can't translate log-linear in my head very well... In the log plots we both agree there seems to be an excess in lower metallicities for the 2-stage formalism (though not the lowest bin); but in my linear plot I don't see it. Am I doing something wrong, or is this just a trick of log vs linear binning?

@ilyamandel @SimonStevenson @jmerritt1

[ilyamandel]

@jeffriley : The log plot is the more useful one in this case, because our metallicity sampling is uniform in the log. If the yield of merging BBHs were independent of metallicity, the output would be flat (up to Poisson fluctuations) on a plot with a log abscissa, but would be dominated by low metallicities on a plot with a linear abscissa. And because many of the low-metallicity logarithmic bins end up being bunched together into a single linear bin, you just see an integrated outcome on the linear abscissa plot, which doesn't allow you to resolve the local features at Z ~ 0.0002.
In any case, I think you, @SimonStevenson , and @jmerritt1 have all confirmed that there is a feature there, and it's still there after my latest changes. Thank you! I will try to understand why.

[ilyamandel]

Heads up: I have made significant changes to the convective envelope mass models in #1230, which will likely impact these tests.

[SimonStevenson]

I reran a test of the merging BBH yield as a function of metallicity using the latest dev (v03.04.00) including #1230 . It seems to have resolved the issue with the over density at low metallicity. Do you have any ideas as to which change in particular was responsible, @ilyamandel ?

Note that these do not include tides, so do not test the implications of those changes for tides. Commands used:

./COMPAS -n 1000000 --metallicity-distribution LOGUNIFORM --common-envelope-formalism ENERGY --initial-mass-min 10.0 --output-container DEFAULT > default.out &
./COMPAS -n 1000000 --metallicity-distribution LOGUNIFORM --common-envelope-formalism TWO_STAGE --initial-mass-min 10.0 --output-container TWO_STAGE > TWO_STAGE.out &

Given this result, I suggest we close this issue.

[ilyamandel]

Hi @SimonStevenson ,

Thank you very much for this test.

The main issue for the 2-stage common envelope was that some stars developed a sudden and deep convective outer envelope at the very start of the HG. This was artificial, unintended behaviour.

Details: The issue was that Tmin (when the convective envelope is supposed to reach 99% of its maximum -- but, in practice, estimated as the COMPAS stellar temperature at the onset of the AGB) is sometimes larger than Tonset (when the envelope just starts becoming convective, estimated according to Eq. 6 of Picker et al.).

For example, for a 5 Msun star at solar metallicity, Tonset = 4144 K (that just depends on metallicity, not mass), but Tmin = 4401 K.

This created a serious problem, because when the current temperature is still much larger than either Tmin or Tonset (say, at the start of the HG), the envelope is already fully convective rather than fully radiative (see Eq. 7 of Picker et al.). Now, looking at Figure 4 of Picker et al., it seems that the trouble is with COMPAS stars being too hot when they start the AGB relative to MESA stars. Not crazy hot, if you think in log_10 (T/K) space. COMPAS says stars start the AGB at log_10 (T/K) ~ 3.64 in this mass range, while Picker et al. figure 4 suggests perhaps log_10 (T/K) ~ 3.54 would be more appropriate according to MESA. This isn't a COMPAS problem -- cf. Figure 1 of Hurley, Pols, Tout, 2000 -- but simply a difference in what the HR diagrams look like between MESA and Pols tracks. But that difference of 0.1 dex in temperature becomes a huge problem for both 2-stage CEs and tides.

The solution I implemented was to adjust T_onset. Just looking at Figure 4 of Picker+, it seems that T_onset is always about 0.1 dex hotter than T_min. So I changed the T_onset to be relative to T_min, and then keeping the rest of the fits the same. It may not be ideal, but it is the simplest fix I can think of, and at least it will always do the sensible thing: T_onset will be guaranteed to be hotter than T_min, and the star will gradually develop a convective envelope as it expands and cools.

I can't tell specifically why it manifested as a peak at low Z, though I suspect it had to do with the fact that the original fit for T_onset in Picker+ was based on metallicity, so at particular metallicities, the relative values of T_min and T_onset changed for particular stellar masses.

However, since the new solution is clearly more sensible and your test, as well as my own tests on individual stars, confirms it, I agree with your suggestion to close this issue.