Reduction in production of BHBH from CHE

[jeffriley]

Significant reduction in production of BHBH from CHE v03.05.02 vs v03.10.01

The following command, using v03.05.02 (*on my server - see notes below):

./compas -n 1 --random-seed 4907 
                       --metallicity 0.0017724805683766772 
                       --initial-mass-1 41.997969179405416 --initial-mass-2 41.997969179405416 
                       --semi-major-axis 0.06926436655793679
                       --chemically-homogeneous-evolution-mode PESSIMISTIC
                       --enhance-CHE-lifetimes-luminosities TRUE
                       --scale-CHE-mass-loss-with-surface-helium-abundance TRUE
                       --enable-rotationally-enhanced-mass-loss TRUE
                       --scale-terminal-wind-velocity-with-metallicity-power 0.0
                       --mass-change-fraction 0.005 --radial-change-fraction 0.005
                       --tides-prescription NONE
                       --pulsational-pair-instability-prescription HENDRIKS
                       --PPI-CO-Core-Shift-Hendriks 1.111111
                       --PPI-upper-limit 80.0 --PISN-lower-limit 80.0
                       --common-envelope-formalism TWO_STAGE
                       --wolf-rayet-multiplier 1.291550 

results in the following output:

COMPAS v03.05.02
Compact Object Mergers: Population Astrophysics and Statistics
by Team COMPAS (http://compas.science/index.html)
A binary star simulator

Go to https://compas.readthedocs.io/en/latest/index.html for the online documentation
Check https://compas.readthedocs.io/en/latest/pages/whats-new.html to see what's new in the latest release

Start generating binaries at Thu Dec  5 09:42:21 2024

0: DCO formed: (Chemically_Homogeneous -> Black_Hole) + (Chemically_Homogeneous -> Black_Hole)

Generated 1 of 1 binaries requested

Simulation completed

End generating binaries at Thu Dec  5 09:42:22 2024

Clock time = 0.101016 CPU seconds
Wall time  = 0000:00:00 (hhhh:mm:ss)

whereas the same command using COMPAS v03.10.01 (again, on my server - see notes below) results in:

COMPAS v03.10.01
Compact Object Mergers: Population Astrophysics and Statistics
by Team COMPAS (http://compas.science/index.html)
A binary star simulator

Go to https://compas.readthedocs.io/en/latest/index.html for the online documentation
Check https://compas.readthedocs.io/en/latest/pages/whats-new.html to see what's new in the latest release

Start generating binaries at Thu Dec  5 09:50:36 2024

0: Stars merged: (Main_Sequence_>_0.7 -> Main_Sequence_>_0.7) + (Main_Sequence_>_0.7 -> Main_Sequence_>_0.7)

Generated 1 of 1 binaries requested

Simulation completed

End generating binaries at Thu Dec  5 09:50:36 2024

Clock time = 0.004449 CPU seconds
Wall time  = 0000:00:00 (hhhh:mm:ss)

In a run of 1 million binaries using the command shown above, on v03.05.02 we had ~650 CHE BHBHs, but on v03.10.01 we had none (see slack). So it would seem we are treating CHE binaries differently in later versions of COMPAS (post v03.05.02), and perhaps more differently than we should...

We should check that the CHE treatment in the latest version of COMPAS is correct.

*on my server (notes above).

If I run the same commands on my desktop, rather than my server, I get different results.

For v03.05.02 I get

COMPAS v03.05.02
Compact Object Mergers: Population Astrophysics and Statistics
by Team COMPAS (http://compas.science/index.html)
A binary star simulator

Go to https://compas.readthedocs.io/en/latest/index.html for the online documentation
Check https://compas.readthedocs.io/en/latest/pages/whats-new.html to see what's new in the latest release

Start generating binaries at Thu Dec  5 10:20:01 2024

0: Unbound binary: (Chemically_Homogeneous -> Neutron_Star) + (Chemically_Homogeneous -> Neutron_Star)

Generated 1 of 1 binaries requested

Simulation completed

End generating binaries at Thu Dec  5 10:20:01 2024

Clock time = 0.030105 CPU seconds
Wall time  = 0000:00:00 (hhhh:mm:ss)

for v03.10.01 I get the same results as the run on my server.

I've noticed this before, and put it down to differences in the compiler implementations, gsl library, os, etc. - and I still think that's where the difference is. For completeness, my server versions are:

OS: Debian GNU/Linux 11 (bullseye)
g++: Debian 10.2.1-6
gsl library (COMPAS uses the gsl random number generator): 2.6+dfsg-2

On my desktop:

OS: Ubuntu 22.04.2 LTS
g++: Ubuntu 11.4.0-1ubuntu1~22.04
gsl library: 2.7.1+dfsg-3

We probably should document this somewhere, and caution users against running on different systems for the same project and combining the results - doing that could result in inconsistencies in the results. It probably only matters for really big runs - I imagine smaller runs would all be done on the same machine. For really big runs that might be distributed across different machines, using the same docker image is probably the safest option - unless the user can be sure each system is identical wrt OS, g++, and gsl library.

[dumbscholar]

With the same command (single ZAMS binary), with version 3.07.02, I get this on M3 Macbook Pro (Sonoma 14.6.1)

g++: Apple clang version 16.0.0 (clang-1600.0.26.4)
gsl: version 2.8

COMPAS v03.07.02
Compact Object Mergers: Population Astrophysics and Statistics
by Team COMPAS (http://compas.science/index.html)
A binary star simulator

Go to https://compas.readthedocs.io/en/latest/index.html for the online documentation
Check https://compas.readthedocs.io/en/latest/pages/whats-new.html to see what's new in the latest release

Start generating binaries at Thu Dec 5 02:01:47 2024

0: DCO formed: (Chemically_Homogeneous -> Black_Hole) + (Chemically_Homogeneous -> Black_Hole)

Generated 1 of 1 binaries requested

Simulation completed

End generating binaries at Thu Dec 5 02:01:47 2024

Clock time = 0.016476 CPU seconds
Wall time = 0000:00:00 (hhhh:mm:ss)

[jeffriley]

@ilyamandel changed behaviour first appears in v03.08.01 (it is not evident in v03.08.00)

[veome22]

@ilyamandel, @jeffriley I may have found the problem. I think the issue was introduced in PR #1279. When checking if a star should be promoted to Type CHE in BaseBinaryStar.cpp, we used to set the stellar rotation to the orbital period first, and then check if the rotation exceeds the CHE rotation threshold. But the newer version directly checks for the spin condition WITHOUT first assuming that the star has spun up to match the orbital frequency. If there is no mechanism to spin up the star (e.g. tides), then the star will never be able to jump to CHE.

(BaseBinaryStar.cpp, Line 293, v3.07.05 -> v3.08.00)

if (OPTIONS->CHEMode() != CHE_MODE::NONE) {                                                                                                         // CHE enabled?

        // CHE enabled, update rotational frequency for constituent stars - assume tidally locked

        double omega = OrbitalAngularVelocity();                                                                                                        // orbital angular velocity

-       m_Star1->SetOmega(omega);
-       m_Star2->SetOmega(omega);

        // check for CHE
        //
        // because we've changed the rotational frequency of the constituent stars we
        // have to reset the stellar type - at this stage, based on their rotational
        // frequency at birth, they may have already been assigned one of MS_LTE_07,
        // MS_GT_07, or CHEMICALLY_HOMOGENEOUS
        //
        // here we need to change from MS_* -> CH, or from CH->MS* based on the
        // newly-assigned rotational frequencies
        // star 1
        if (utils::Compare(m_Star1->Omega(), m_Star1->OmegaCHE()) >= 0) {                                                                               // star 1 CH?
            if (m_Star1->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star1->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);    // yes, switch if not already Chemically Homogeneous
+           m_Star1->SetOmega(omega);
        }
        else if (m_Star1->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
            if (m_Star1->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.7 Msol - switch if necessary
        }
        else {
            if (m_Star1->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star1->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
        }
        // star 2
        if (utils::Compare(m_Star2->Omega(), m_Star2->OmegaCHE()) >= 0) {                                                                               // star 2 CH?
            if (m_Star2->StellarType() != STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) (void)m_Star2->SwitchTo(STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS, true);    // yes, switch if not already Chemically Homogeneous
+           m_Star2->SetOmega(omega);
        }
        else if (m_Star2->MZAMS() <= 0.7) {                                                                                                             // no - MS - initial mass determines actual type  (don't use utils::Compare() here)
            if (m_Star2->StellarType() != STELLAR_TYPE::MS_LTE_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_LTE_07, true);                              // MS <= 0.0 Msol - switch if necessary
        }
        else {
            if (m_Star2->StellarType() != STELLAR_TYPE::MS_GT_07) (void)m_Star2->SwitchTo(STELLAR_TYPE::MS_GT_07, true);                                // MS > 0.7 Msol - switch if necessary
        }

I think the intended solution was to let ProcessTides() spin up the stars, but that first requires the stars to be of type CHE.

ProcessTides():

case TIDES_PRESCRIPTION::NONE: {                                                                                    // NONE - tides not enabled
            
                // do nothing, except for CHE stars which are allowed to remain CHE even though this does not preserve angular momentum
                if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega);
                if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega);
                    
            } break;

The previous approach of always setting stellar rotation rates to the orbital frequency was obviously problematic, because it interfered with spins from tides. But we need some other way of deciding to switch to CHE. @ilyamandel, you had suggested using an orbital period threshold, such that all ZAMS binaries below, say 3 days, would automatically be considered chemically homogeneously evolving, and thus in synchronous and circular orbits. We can do this as an override, so that it applies no matter what tidal prescription was chosen.