Skip to content

Commit 5856d88

Browse files
authored
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
2 parents 4f8a2f4 + ffb4250 commit 5856d88

File tree

17 files changed

+210
-118
lines changed

17 files changed

+210
-118
lines changed

compas_matlab_utils/ComparisonPlots.m

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -258,8 +258,9 @@ function HMXBplot(file, name, fignumberHMXB, colour, point)
258258
roche2Switch=h5read(file,'/BSE_Switch_Log/RocheLobe(2)');
259259
radius1Switch=h5read(file,'/BSE_Switch_Log/Radius(1)');
260260
radius2Switch=h5read(file,'/BSE_Switch_Log/Radius(2)');
261-
relevantBinaryFrom1 = whichSwitch==1 & fromSwitch==1 & toSwitch==2 & star2Switch==14 & radius1Switch>0.8*roche1Switch & M1Switch>15;
262-
relevantBinaryFrom2 = whichSwitch==2 & fromSwitch==1 & toSwitch==2 & star1Switch==14 & radius2Switch>0.8*roche2Switch & M2Switch>15;
261+
isMergerSwitch=h5read(file,'/BSE_Switch_Log/Is_Merger');
262+
relevantBinaryFrom1 = whichSwitch==1 & fromSwitch==1 & toSwitch==2 & star2Switch==14 & radius1Switch>0.8*roche1Switch & M1Switch>15 & ~isMergerSwitch;
263+
relevantBinaryFrom2 = whichSwitch==2 & fromSwitch==1 & toSwitch==2 & star1Switch==14 & radius2Switch>0.8*roche2Switch & M2Switch>15 & ~isMergerSwitch;
263264
MBH=[MBH; M2Switch(relevantBinaryFrom1); M1Switch(relevantBinaryFrom2)];
264265
MO=[MO; M1Switch(relevantBinaryFrom2); M2Switch(relevantBinaryFrom1)];
265266
figure(fignumberHMXB), scatter(MBH, MO, point, 'filled', colour, 'DisplayName', name); hold on;
@@ -342,7 +343,8 @@ function DWDplot(file, name, fignumberDWD, colourNMT, colourMT, point)
342343
M2Switch=h5read(file, '/BSE_Switch_Log/Mass(2)');
343344
aSwitch=h5read(file, '/BSE_Switch_Log/SemiMajorAxis');
344345
SeedSwitch=h5read(file, '/BSE_Switch_Log/SEED');
345-
WDIndex=find(Type1>=10 & Type1<=12 & Type2>=10 & Type2<=12);
346+
isMergerSwitch=h5read(file,'/BSE_Switch_Log/Is_Merger');
347+
WDIndex=find(Type1>=10 & Type1<=12 & Type2>=10 & Type2<=12 & ~isMergerSwitch);
346348
ind=unique(SeedSwitch(WDIndex));
347349
[~,ZAMSIndex]=ismember(ind,SeedZAMS);
348350
MAZAMS=(M1ZAMS(ZAMSIndex)+M2ZAMS(ZAMSIndex)).*aZAMS(ZAMSIndex);

compas_matlab_utils/CosmicHistoryIntegrator.m

Lines changed: 31 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,13 @@
1-
function [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection, x]=...
1+
function [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
2+
MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
23
CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
34
% Integrator for the binary black hole merger rate over cosmic history
45
% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics)
56
% software package
67
%
78
% USAGE:
8-
% [Zlist, MergerRateByRedshiftByZ, SFR, Zweight, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
9+
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
10+
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]]=...
911
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated [,makeplots])
1012
%
1113
% INPUTS:
@@ -19,32 +21,44 @@
1921
% makeplots: if set to 1, generates a set of useful plots (default = 0)
2022
%
2123
% OUTPUTS:
24+
% SFR is a vector of size length(zlistformation) containing the star formation rate
25+
% (solar masses per Mpc^3 of comoving volume per year of source time)
2226
% Zlist is a vector of metallicities, taken from the COMPAS run input file
27+
% Mtlist is a list of total mass bins
28+
% etalist is a list of symmetric mass ratio bins
29+
% FormationRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
30+
% which contains a formation rate of merging double compact objects in the given redshift
31+
% and metallicity bin, in units of formed DCOs per Mpc^3 of comoving volume per
32+
% year of source time
33+
% FormationRateByRedshiftByMtByEta is a matrix of size length(zformationlist)
34+
% X length(Mtlist) X length(etalist) which contains a formation rate of merging double compact objects
35+
% in the given redshift, total mass and eta bin, in units of formed DCOs per Mpc^3
36+
% of comoving volume per year of source time
2337
% MergerRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
24-
% which contains a merger rate of merging compact objects in the given redshift
38+
% which contains a merger rate of double compact objects in the given redshift
2539
% and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
2640
% year of source time
27-
% MergerRateByRedshiftByMtByEta is a matrix of size llength(zformationlist)
28-
% X length(Mtlist) X length(etalist) which contains a merger rate of merging compact objects
41+
% MergerRateByRedshiftByMtByEta is a matrix of size length(zformationlist)
42+
% X length(Mtlist) X length(etalist) which contains a merger rate of double compact objects
2943
% in the given redshift, total mass and eta bin, in units of mergers per Mpc^3
3044
% of comoving volume per year of source time
31-
% Mtlist is a list of total mass bins
32-
% etalist is a list of symmetric mass ratio bins
33-
% SFR is a vector of size length(zlistformation) containing the star formation rate
34-
% (solar masses per Mpc^3 of comoving volume per year of source time)
45+
% zlistdetection is a vector of redshifts at which detection rates are
46+
% computed (a subset of zlistformation going up to zmaxdetection)
3547
% Rdetection is a matrix of size length(zlistdetection) X length(Mtlist) X
3648
% length(etalist) containing the detection rate per year of observer time
3749
% from a given redshift bin and total mass and symmetric mass ratio pixel
3850
% DetectableMergerRate is a matrix of the same size as Rdetection but
3951
% containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
4052
% volume per year of source time
41-
% zlistdetection is a vector of redshifts at which detection rates are
42-
% computed (a subset of zlistformation going up to zmaxdetection)
53+
4354
%
4455
% EXAMPLE:
4556
% zlist=0:0.01:10;
46-
% [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection,] = ...
57+
% [SFR, Zlist, Mtlist, etalist, FormationRateByRedshiftByZ, FormationRateByRedshiftByMtByEta, ...
58+
% MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, zlistdetection, Rdetections, DetectableMergerRate]=...
4759
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
60+
% figure(10), semilogy(zlist, sum(MergerRateByRedshiftByZ,2)*1e9,'LineWidth',3), set(gca,'FontSize',20),
61+
% xlabel('Redshift z'), ylabel('Formation rate of merging DCO per Gpc^3 per yr')
4862
%
4963

5064

@@ -76,19 +90,23 @@
7690
dz=zlistformation(2)-zlistformation(1);
7791
etalist=0.01:0.01:0.25;
7892
Mtlist=1:1:ceil(max(M1+M2));
93+
FormationRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
94+
FormationRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
7995
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
8096
MergerRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
8197
x=zeros(size(M1));
8298
for(i=1:length(M1)),
8399
Zcounter=find(Zlist>=Z(i),1);
84100
eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
85101
etaindex=ceil(eta*100);
102+
Mtindex=ceil(M1(i)+M2(i));
103+
FormationRateByRedshiftByZ(:,Zcounter)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
104+
FormationRateByRedshiftByMtByEta(:,Mtindex,etaindex)=transpose(SFR).*Zweight(:,Zcounter)/Msimulated;
86105
tLform=tL+Tdelay(i); %lookback time of when binary would have to form in order to merge at lookback time tL
87106
firsttooearlyindex=find((tLform)>max(tL),1);
88107
if(isempty(firsttooearlyindex)), firsttooearlyindex=length(tL)+1; end;
89108
zForm=interp1(tL,zlistformation,tLform(1:firsttooearlyindex-1));
90109
zFormindex=ceil((zForm-zlistformation(1))./dz)+1;
91-
Mtindex=ceil(M1(i)+M2(i));
92110
if(~isempty(zFormindex))
93111
x(i)=SFR(zFormindex(1))*Zweight(zFormindex(1),Zcounter)/Msimulated;
94112
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)=...
@@ -98,17 +116,6 @@
98116
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) + ...
99117
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
100118
end;
101-
%for(k=1:length(zlistformation)), %merger redshift index
102-
% if((Tdelay(i)+tL(k)) > max(tL)), continue; end; %binary can't merge this early
103-
% zformindex=find(tL>=(Tdelay(i)+tL(k)),1);
104-
% Mtzindex=ceil((M1(i)+M2(i))*(1+zlistformation(k)));
105-
% MergerRateByRedshiftByZ(k,Zcounter)=...
106-
% MergerRateByRedshiftByZ(k,Zcounter)+...
107-
% SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
108-
% MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) =...
109-
% MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) + ...
110-
% SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
111-
%end;
112119
end;
113120

114121
zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));

compas_python_utils/preprocessing/compasConfigDefault.yaml

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
##~!!~## COMPAS option values
2-
##~!!~## File Created Thu Jul 17 09:38:18 2025 by COMPAS v03.21.00
2+
##~!!~## File Created Fri Jul 18 08:19:00 2025 by COMPAS v03.21.00
33
##~!!~##
44
##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
55
##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -27,10 +27,10 @@ booleanChoices:
2727
# --check-photon-tiring-limit: False # Default: False
2828
# --use-mass-loss: True # Default: True
2929
# --enable-rotationally-enhanced-mass-loss: False # Default: False
30-
# --enhance-CHE-lifetimes-luminosities: False # Default: False
30+
# --enhance-CHE-lifetimes-luminosities: True # Default: True
3131
# --expel-convective-envelope-above-luminosity-threshold: False # Default: False
3232
# --natal-kick-for-PPISN: False # Default: False
33-
# --scale-CHE-mass-loss-with-surface-helium-abundance: False # Default: False
33+
# --scale-CHE-mass-loss-with-surface-helium-abundance: True # Default: True
3434

3535
### BINARY PROPERTIES
3636
# --allow-touching-at-birth: False # Default: False # record binaries that have stars touching at birth in output files
@@ -166,6 +166,7 @@ numericalChoices:
166166
# --common-envelope-mass-accretion-max: 0.100000 # Default: 0.100000 # For 'MACLEOD+2014' [Msol]
167167
# --common-envelope-mass-accretion-min: 0.040000 # Default: 0.040000 # For 'MACLEOD+2014' [Msol]
168168
# --common-envelope-recombination-energy-density: 15000000000000.00 # Default: 15000000000000.00
169+
# --common-envelope-second-stage-beta: 0.000000 # Default: 0.000000
169170
# --common-envelope-slope-kruckow: -0.833333 # Default: -0.833333
170171
# --maximum-mass-donor-nandez-ivanova: 2.000000 # Default: 2.000000
171172

@@ -197,7 +198,8 @@ numericalChoices:
197198
# --mcbur1: 1.600000 # Default: 1.600000
198199
# --muller-mandel-kick-multiplier-BH: 200.00 # Default: 200.00 # scaling prefactor for BH kicks when using the 'MULLERMANDEL' kick magnitude distribution
199200
# --muller-mandel-kick-multiplier-NS: 520.00 # Default: 520.00 # scaling prefactor for NS kicks when using the 'MULLERMANDEL' kick magnitude distribution
200-
# --muller-mandel-sigma-kick: 0.300000 # Default: 0.300000 # kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
201+
# --muller-mandel-sigma-kick-BH: 0.300000 # Default: 0.300000 # BH kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
202+
# --muller-mandel-sigma-kick-NS: 0.300000 # Default: 0.300000 # NS kick scatter when using the 'MULLERMANDEL' kick magnitude distribution
201203
# --neutrino-mass-loss-BH-formation-value: 0.100000 # Default: 0.100000
202204
# --pisn-lower-limit: 60.000000 # Default: 60.000000 # Minimum core mass for PISN [Msol]
203205
# --pisn-upper-limit: 135.00 # Default: 135.00 # Maximum core mass for PISN [Msol]
@@ -272,6 +274,7 @@ stringChoices:
272274
# --common-envelope-formalism: 'ENERGY' # Default: 'ENERGY' # Options: ['TWO_STAGE','ENERGY']
273275
# --common-envelope-lambda-prescription: 'LAMBDA_NANJING' # Default: 'LAMBDA_NANJING' # Options: ['LAMBDA_DEWI','LAMBDA_KRUCKOW','LAMBDA_NANJING','LAMBDA_LOVERIDGE','LAMBDA_FIXED'] # Xu & Li 2010
274276
# --common-envelope-mass-accretion-prescription: 'ZERO' # Default: 'ZERO' # Options: ['CHEVALIER','MACLEOD','UNIFORM','CONSTANT','ZERO']
277+
# --common-envelope-second-stage-gamma-prescription: 'ISOTROPIC' # Default: 'ISOTROPIC' # Options: ['ARBITRARY','MACLEOD_LINEAR','CIRCUMBINARY','ISOTROPIC','JEANS']
275278

276279
### TIDES
277280
# --tides-prescription: 'NONE' # Default: 'NONE' # Options: ['KAPIL2025','PERFECT','NONE']

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

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -398,7 +398,7 @@ Default = FALSE
398398

399399
**--enhance-CHE-lifetimes-luminosities** |br|
400400
Enhance lifetimes and luminosities of CH stars using a fit to detailed models from Szecsi et al. (2015)
401-
Default = FALSE
401+
Default = TRUE
402402

403403
**--envelope-state-prescription** |br|
404404
Prescription for determining whether the envelope of the star is convective or radiative. |br|
@@ -513,7 +513,7 @@ Default = Sampled from IMF
513513

514514
**--initial-mass-2** |br|
515515
Initial mass for the secondary star when evolving in BSE mode (:math:`M_\odot`). |br|
516-
Default = Sampled from IMF
516+
Default = Sampled from the mass ratio distribution specified by ``--mass-ratio-distribution`` (see also ``--mass-ratio-max``, ``-mass-ratio-min``, ``--minimum-secondary-mass'')
517517
518518
**--initial-mass-function [ -i ]** |br|
519519
Initial mass function. |br|
@@ -967,8 +967,12 @@ Default = 200.0
967967
Scaling prefactor for NS kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
968968
Default = 520.0
969969

970-
**--muller-mandel-sigma-kick** |br|
971-
Scatter width for NS and BH kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
970+
**--muller-mandel-sigma-kick-BH** |br|
971+
Scatter width for BH kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
972+
Default = 0.3
973+
974+
**--muller-mandel-sigma-kick-NS** |br|
975+
Scatter width for NS kicks when using the `MULLERMANDEL` kick magnitude distribution |br|
972976
Default = 0.3
973977

974978
.. _options-props-N:
@@ -1290,7 +1294,7 @@ Default = DECIN2023 |br|
12901294
**--scale-CHE-mass-loss-with-surface-helium-abundance** |br|
12911295
Scale mass loss for chemically homogeneously evolving (CHE) stars with the surface helium abundance.
12921296
Transition from OB to WR mass loss towards the end of the main sequence.
1293-
Default = False
1297+
Default = TRUE
12941298

12951299
**--scale-terminal-wind-velocity-with-metallicity-power** |br|
12961300
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
15641568
--kick-magnitude-distribution, --kick-magnitude-sigma-CCSN-BH, --kick-magnitude-sigma-CCSN-NS, --kick-magnitude-sigma-ECSN, --kick-magnitude-sigma-USSN,
15651569
--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,
15661570
--kick-magnitude-random, --kick-magnitude-random-1, --kick-magnitude-random-2, --kick-scaling-factor, -muller-mandel-kick-multiplier-BH,
1567-
--muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick
1568-
1569-
--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
1571+
--muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick-BH, --muller-mandel-sigma-kick-NS, --kick-direction,
1572+
--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
15701573

15711574
:ref:`Back to Top <options-props-top>`
15721575

online-docs/pages/whats-new.rst

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,13 @@ What's new
33

44
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``.
55

6-
**03.20=1.00 July 17, 2025**
6+
**03.22.00 July 18, 2025**
7+
* Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
8+
* 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``)
9+
* 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)
10+
* Fixed the calculation of Mandel & Muller kicks; split ``--muller-mandel-sigma-kick`` into ``--muller-mandel-sigma-kick-NS`` and ``--muller-mandel-sigma-kick-BH``
11+
12+
**03.21.00 July 17, 2025**
713

814
* Deprecated mass loss prescription MERRITT2024 in favour of MERRITT2025
915
* Added version strings for gsl, boost, and HDF5 to COMPAS splashscreen

0 commit comments

Comments
 (0)