Merge pull request #1409 from TeamCOMPAS/Matlab

03.22.00: fixes to nuclear timescale mass transfer, new 2-stage CE options, etc:

Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
Added options to set beta and gamma prescription for second stage of 2-stage CE (--common-envelope-second-stage-beta, --common-envelope-second-stage-gamma-prescription)
Fixed a bug in CalculateZetaEquilibrium(), which impacted when mass transfer is declared nuclear (and how conservative it is)
Now calculate mass accretion rate for nuclear timescale mass transfer on the fly to match with donor mass loss rate set by donor mass loss (required to fit into Roche lobe) divided by time step
Fixed random draws of SN kicks to avoid artificial pile-up at boundaries of distribution
Split --muller-mandel-sigma-kick into --muller-mandel-sigma-kick-NS and --muller-mandel-sigma-kick-BH

Author: ilyamandel

compas_matlab_utils/ComparisonPlots.m

@@ -258,8 +258,9 @@ function HMXBplot(file, name, fignumberHMXB, colour, point)
     roche2Switch=h5read(file,'/BSE_Switch_Log/RocheLobe(2)');
     radius1Switch=h5read(file,'/BSE_Switch_Log/Radius(1)');
     radius2Switch=h5read(file,'/BSE_Switch_Log/Radius(2)');
-    relevantBinaryFrom1 = whichSwitch==1 & fromSwitch==1 & toSwitch==2 & star2Switch==14 & radius1Switch>0.8*roche1Switch & M1Switch>15;
-    relevantBinaryFrom2 = whichSwitch==2 & fromSwitch==1 & toSwitch==2 & star1Switch==14 & radius2Switch>0.8*roche2Switch & M2Switch>15;
+    isMergerSwitch=h5read(file,'/BSE_Switch_Log/Is_Merger');
+    relevantBinaryFrom1 = whichSwitch==1 & fromSwitch==1 & toSwitch==2 & star2Switch==14 & radius1Switch>0.8*roche1Switch & M1Switch>15 & ~isMergerSwitch;
+    relevantBinaryFrom2 = whichSwitch==2 & fromSwitch==1 & toSwitch==2 & star1Switch==14 & radius2Switch>0.8*roche2Switch & M2Switch>15 & ~isMergerSwitch;
     MBH=[MBH; M2Switch(relevantBinaryFrom1); M1Switch(relevantBinaryFrom2)];
     MO=[MO; M1Switch(relevantBinaryFrom2); M2Switch(relevantBinaryFrom1)];
     figure(fignumberHMXB), scatter(MBH, MO, point, 'filled', colour, 'DisplayName', name); hold on;
@@ -342,7 +343,8 @@ function DWDplot(file, name, fignumberDWD, colourNMT, colourMT, point)
     M2Switch=h5read(file, '/BSE_Switch_Log/Mass(2)');
     aSwitch=h5read(file, '/BSE_Switch_Log/SemiMajorAxis');
     SeedSwitch=h5read(file, '/BSE_Switch_Log/SEED');
-    WDIndex=find(Type1>=10 & Type1<=12 & Type2>=10 & Type2<=12);
+    isMergerSwitch=h5read(file,'/BSE_Switch_Log/Is_Merger');
+    WDIndex=find(Type1>=10 & Type1<=12 & Type2>=10 & Type2<=12 & ~isMergerSwitch);
     ind=unique(SeedSwitch(WDIndex));
     [~,ZAMSIndex]=ismember(ind,SeedZAMS);
     MAZAMS=(M1ZAMS(ZAMSIndex)+M2ZAMS(ZAMSIndex)).*aZAMS(ZAMSIndex);

compas_matlab_utils/CosmicHistoryIntegrator.m

@@ -1,11 +1,13 @@
-function [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta,  Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection, x]=...
+function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
+    MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
     CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
 % Integrator for the binary black hole merger rate over cosmic history
 % COMPAS (Compact Object Mergers: Population Astrophysics and Statistics) 
 % software package
 %
 % USAGE: 
-% [Zlist, MergerRateByRedshiftByZ, SFR, Zweight, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
+% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
+%    MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]]=...
 %    CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated [,makeplots])
 %
 % INPUTS:
@@ -19,32 +21,44 @@
 %   makeplots:  if set to 1, generates a set of useful plots (default = 0)
 %
 % OUTPUTS: 
+%   SFR is a vector of size length(zlistformation) containing the star formation rate 
+% (solar masses per Mpc^3 of comoving volume per year of source time)
 %   Zlist is a vector of metallicities, taken from the COMPAS run input file
+%   Mtlist is a list of total mass bins
+%   etalist is a list of symmetric mass ratio bins
+%   FormationRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist) 
+% which contains a formation rate of merging double compact objects in the given redshift 
+% and metallicity bin, in units of formed DCOs per Mpc^3 of comoving volume per
+% year of source time
+%   FormationRateByRedshiftByMtByEta is a matrix of size length(zformationlist) 
+% X length(Mtlist) X length(etalist) which contains a formation rate of merging double compact objects 
+% in the given redshift, total mass and eta bin, in units of formed DCOs per Mpc^3 
+% of comoving volume per year of source time
 %   MergerRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist) 
-% which contains a merger rate of merging compact objects in the given redshift 
+% which contains a merger rate of double compact objects in the given redshift 
 % and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
 % year of source time
-%   MergerRateByRedshiftByMtByEta is a matrix of size llength(zformationlist) 
-% X length(Mtlist) X length(etalist) which contains a merger rate of merging compact objects 
+%   MergerRateByRedshiftByMtByEta is a matrix of size length(zformationlist) 
+% X length(Mtlist) X length(etalist) which contains a merger rate of double compact objects 
 % in the given redshift, total mass and eta bin, in units of mergers per Mpc^3 
 % of comoving volume per year of source time
-%   Mtlist is a list of total mass bins
-%   etalist is a list of symmetric mass ratio bins
-%   SFR is a vector of size length(zlistformation) containing the star formation rate 
-% (solar masses per Mpc^3 of comoving volume per year of source time)
+%   zlistdetection is a vector of redshifts at which detection rates are
+% computed (a subset of zlistformation going up to zmaxdetection)
 %   Rdetection is a matrix of size length(zlistdetection) X length(Mtlist) X
 % length(etalist) containing the detection rate per year of observer time
 % from a given redshift bin and total mass and symmetric mass ratio pixel
 %   DetectableMergerRate is a matrix of the same size as Rdetection but 
 % containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
 % volume per year of source time
-%   zlistdetection is a vector of redshifts at which detection rates are
-% computed (a subset of zlistformation going up to zmaxdetection)
+
 %
 % EXAMPLE:
 % zlist=0:0.01:10;
-% [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta,  Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection,] = ...
+% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
+%    MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
 % CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
+% figure(10), semilogy(zlist, sum(MergerRateByRedshiftByZ,2)*1e9,'LineWidth',3), set(gca,'FontSize',20),
+% xlabel('Redshift z'), ylabel('Formation rate of merging DCO per Gpc^3 per yr')
 % 
 
 
@@ -76,19 +90,23 @@
 dz=zlistformation(2)-zlistformation(1);
 etalist=0.01:0.01:0.25;
 Mtlist=1:1:ceil(max(M1+M2));
+FormationRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
+FormationRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
 MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
 MergerRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
 x=zeros(size(M1));
 for(i=1:length(M1)),
     Zcounter=find(Zlist>=Z(i),1);
     eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
     etaindex=ceil(eta*100);
+    Mtindex=ceil(M1(i)+M2(i));
+    FormationRateByRedshiftByZ(:,Zcounter)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
+    FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
     tLform=tL+Tdelay(i);    %lookback time of when binary would have to form in order to merge at lookback time tL
     firsttooearlyindex=find((tLform)>max(tL),1);
     if(isempty(firsttooearlyindex)), firsttooearlyindex=length(tL)+1; end;
     zForm=interp1(tL,zlistformation,tLform(1:firsttooearlyindex-1));
     zFormindex=ceil((zForm-zlistformation(1))./dz)+1;
-    Mtindex=ceil(M1(i)+M2(i));
     if(~isempty(zFormindex))
         x(i)=SFR(zFormindex(1))*Zweight(zFormindex(1),Zcounter)/Msimulated;
         MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)=...
@@ -98,17 +116,6 @@
             MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) + ...
             transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
     end;
-    %for(k=1:length(zlistformation)),    %merger redshift index
-    %    if((Tdelay(i)+tL(k)) > max(tL)), continue; end;    %binary can't merge this early
-    %    zformindex=find(tL>=(Tdelay(i)+tL(k)),1);
-    %    Mtzindex=ceil((M1(i)+M2(i))*(1+zlistformation(k)));
-    %    MergerRateByRedshiftByZ(k,Zcounter)=...
-    %        MergerRateByRedshiftByZ(k,Zcounter)+...
-    %        SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated; 
-    %    MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) =...
-    %        MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) + ...
-    %        SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
-    %end;
 end;
 
 zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Thu Jul 17 09:38:18 2025 by COMPAS v03.21.00
+##~!!~## File Created Fri Jul 18 08:19:00 2025 by COMPAS v03.21.00
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -27,10 +27,10 @@ booleanChoices:
 #    --check-photon-tiring-limit: False                                    # Default: False
 #    --use-mass-loss: True                                                 # Default: True
 #    --enable-rotationally-enhanced-mass-loss: False                       # Default: False
-#    --enhance-CHE-lifetimes-luminosities: False                           # Default: False
+#    --enhance-CHE-lifetimes-luminosities: True                            # Default: True
 #    --expel-convective-envelope-above-luminosity-threshold: False         # Default: False
 #    --natal-kick-for-PPISN: False                                         # Default: False
-#    --scale-CHE-mass-loss-with-surface-helium-abundance: False            # Default: False
+#    --scale-CHE-mass-loss-with-surface-helium-abundance: True             # Default: True
 
     ### BINARY PROPERTIES
 #    --allow-touching-at-birth: False                                      # Default: False                                                                                                                                                                            # record binaries that have stars touching at birth in output files
@@ -166,6 +166,7 @@ numericalChoices:
 #    --common-envelope-mass-accretion-max: 0.100000                        # Default: 0.100000                                                                                                                                                                         # For 'MACLEOD+2014' [Msol]
 #    --common-envelope-mass-accretion-min: 0.040000                        # Default: 0.040000                                                                                                                                                                         # For 'MACLEOD+2014' [Msol]
 #    --common-envelope-recombination-energy-density: 15000000000000.00     # Default: 15000000000000.00
+#    --common-envelope-second-stage-beta: 0.000000                         # Default: 0.000000
 #    --common-envelope-slope-kruckow: -0.833333                            # Default: -0.833333
 #    --maximum-mass-donor-nandez-ivanova: 2.000000                         # Default: 2.000000
 
@@ -197,7 +198,8 @@ numericalChoices:
 #    --mcbur1: 1.600000                                                    # Default: 1.600000
 #    --muller-mandel-kick-multiplier-BH: 200.00                            # Default: 200.00                                                                                                                                                                           # scaling prefactor for BH kicks when using the 'MULLERMANDEL' kick magnitude distribution
 #    --muller-mandel-kick-multiplier-NS: 520.00                            # Default: 520.00                                                                                                                                                                           # scaling prefactor for NS kicks when using the 'MULLERMANDEL' kick magnitude distribution
-#    --muller-mandel-sigma-kick: 0.300000                                  # Default: 0.300000                                                                                                                                                                         # kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
+#    --muller-mandel-sigma-kick-BH: 0.300000                               # Default: 0.300000                                                                                                                                                                         # BH kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
+#    --muller-mandel-sigma-kick-NS: 0.300000                               # Default: 0.300000                                                                                                                                                                         # NS kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
 #    --neutrino-mass-loss-BH-formation-value: 0.100000                     # Default: 0.100000
 #    --pisn-lower-limit: 60.000000                                         # Default: 60.000000                                                                                                                                                                        # Minimum core mass for PISN [Msol]
 #    --pisn-upper-limit: 135.00                                            # Default: 135.00                                                                                                                                                                           # Maximum core mass for PISN [Msol]
@@ -272,6 +274,7 @@ stringChoices:
 #    --common-envelope-formalism: 'ENERGY'                                 # Default: 'ENERGY'                      # Options: ['TWO_STAGE','ENERGY']
 #    --common-envelope-lambda-prescription: 'LAMBDA_NANJING'               # Default: 'LAMBDA_NANJING'              # Options: ['LAMBDA_DEWI','LAMBDA_KRUCKOW','LAMBDA_NANJING','LAMBDA_LOVERIDGE','LAMBDA_FIXED']                                                     # Xu & Li 2010
 #    --common-envelope-mass-accretion-prescription: 'ZERO'                 # Default: 'ZERO'                        # Options: ['CHEVALIER','MACLEOD','UNIFORM','CONSTANT','ZERO']
+#    --common-envelope-second-stage-gamma-prescription: 'ISOTROPIC'        # Default: 'ISOTROPIC'                   # Options: ['ARBITRARY','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS']
 
     ### TIDES
 #    --tides-prescription: 'NONE'                                          # Default: 'NONE'                        # Options: ['KAPIL2025','PERFECT','NONE']

online-docs/pages/User guide/Program options/program-options-list-defaults.rst

@@ -398,7 +398,7 @@ Default = FALSE
 
 **--enhance-CHE-lifetimes-luminosities** |br|
 Enhance lifetimes and luminosities of CH stars using a fit to detailed models from Szecsi et al. (2015)
-Default = FALSE
+Default = TRUE
 
 **--envelope-state-prescription** |br|
 Prescription for determining whether the envelope of the star is convective or radiative. |br|
@@ -513,7 +513,7 @@ Default = Sampled from IMF
 
 **--initial-mass-2** |br|
 Initial mass for the secondary star when evolving in BSE mode (:math:`M_\odot`). |br|
-Default = Sampled from IMF
+Default = Sampled from the mass ratio distribution specified by ``--mass-ratio-distribution`` (see also ``--mass-ratio-max``, ``-mass-ratio-min``, ``--minimum-secondary-mass'')
 
 **--initial-mass-function [ -i ]** |br|
 Initial mass function. |br|
@@ -967,8 +967,12 @@ Default = 200.0
 Scaling prefactor for NS kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
 Default = 520.0
 
-**--muller-mandel-sigma-kick** |br|
-Scatter width for NS and BH kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
+**--muller-mandel-sigma-kick-BH** |br|
+Scatter width for BH kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
+Default = 0.3
+
+**--muller-mandel-sigma-kick-NS** |br|
+Scatter width for NS kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
 Default = 0.3
 
 .. _options-props-N:
@@ -1290,7 +1294,7 @@ Default = DECIN2023 |br|
 **--scale-CHE-mass-loss-with-surface-helium-abundance** |br|
 Scale mass loss for chemically homogeneously evolving (CHE) stars with the surface helium abundance. 
 Transition from OB to WR mass loss towards the end of the main sequence.
-Default = False
+Default = TRUE
 
 **--scale-terminal-wind-velocity-with-metallicity-power** |br|
 Scale terminal wind velocity with metallicity to this power
@@ -1564,9 +1568,8 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 --kick-magnitude-distribution, --kick-magnitude-sigma-CCSN-BH, --kick-magnitude-sigma-CCSN-NS, --kick-magnitude-sigma-ECSN, --kick-magnitude-sigma-USSN, 
 --black-hole-kicks, --black-hole-kicks-mode, --fix-dimensionless-kick-magnitude, --kick-magnitude, --kick-magnitude-1, --kick-magnitude-2, --kick-magnitude-min, --kick-magnitude-max, 
 --kick-magnitude-random, --kick-magnitude-random-1, --kick-magnitude-random-2, --kick-scaling-factor, -muller-mandel-kick-multiplier-BH, 
---muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick
-
---kick-direction, --kick-direction-distribution, --kick-direction-power, --kick-mean-anomaly-1, --kick-mean-anomaly-2, --kick-phi-1, --kick-phi-2, --kick-theta-1, --kick-theta-2
+--muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick-BH, --muller-mandel-sigma-kick-NS, --kick-direction, 
+--kick-direction-distribution, --kick-direction-power, --kick-mean-anomaly-1, --kick-mean-anomaly-2, --kick-phi-1, --kick-phi-2, --kick-theta-1, --kick-theta-2
 
 :ref:`Back to Top <options-props-top>`
 

online-docs/pages/whats-new.rst

@@ -3,7 +3,13 @@ What's new
 
 Following is a brief list of important updates to the COMPAS code.  A complete record of changes can be found in the file ``changelog.h``.
 
-**03.20=1.00 July 17, 2025**
+**03.22.00 July 18, 2025**
+* Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
+* Added options to set beta and gamma prescription for second stage of 2-stage CE (``--common-envelope-second-stage-beta``, ``--common-envelope-second-stage-gamma-prescription``)
+* Fixed a bug in the calculation of zeta_equilibrium, which impacts when mass transfer is declared to proceed on a nuclear timescale (and hence how conservative it is)
+* Fixed the calculation of Mandel & Muller kicks; split ``--muller-mandel-sigma-kick`` into ``--muller-mandel-sigma-kick-NS`` and ``--muller-mandel-sigma-kick-BH``
+
+**03.21.00 July 17, 2025**
 
 * Deprecated mass loss prescription MERRITT2024 in favour of MERRITT2025
 * Added version strings for gsl, boost, and HDF5 to COMPAS splashscreen