|
| 1 | +function [Zlist, MergerRateByRedshiftByZ, SFRfractionZ]=... |
| 2 | + CosmicHistoryIntegrator(filename,zlist,metallicitychoice, makeplots) |
| 3 | +% Integrator for the binary black hole merger rate over cosmic history |
| 4 | +% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics) |
| 5 | +% software package |
| 6 | +% |
| 7 | +% USAGE: |
| 8 | +% [Zlist, MergerRateByRedshiftByZ]=... |
| 9 | +% CosmicHistoryIntegrator(filename, zlistformation, zlistmerger, zlistdetection, Msimulated, [,makeplots]) |
| 10 | +% |
| 11 | +% INPUTS: |
| 12 | +% filename: name of population synthesis input file |
| 13 | +% should be in COMPAS output h5 format |
| 14 | +% zlistformation: vector of redshifts at which the formation rate is |
| 15 | +% computed |
| 16 | +% zlistmerger: vector of redshifts at which the merger rate is computed |
| 17 | +% zlistdetection: vector of redshifts at which the detection rate is |
| 18 | +% computed |
| 19 | +% Msimulated: total star forming mass represented by the simulation (for |
| 20 | +% normalisation) |
| 21 | +% makeplots: if set to 1, generates a set of useful plots (default = 0) |
| 22 | +% |
| 23 | +% Zlist is a vector of metallicities, taken from input file |
| 24 | +% MergerRateByRedshiftByZ is a matrix of size length(Zlist) X length(zlist) |
| 25 | +% which contains a merger rate of binary black holes in the given redshift |
| 26 | +% and metallicity bin, in units of mergers per Gpc^3 of comoving volume per |
| 27 | +% year of source time |
| 28 | +% |
| 29 | +% EXAMPLE: |
| 30 | +% zlist=0:0.1:5; |
| 31 | +% filename=['~/Work/COMPASresults/popsynth/runs/',... |
| 32 | +% '20170628-Coen-Pessimistic/AllMergers.dat']; |
| 33 | +% [Zlist,MergerRateByRedshiftByZ]=CosmicHistoryIntegrator(filename,zlist); |
| 34 | +% figure(1),colormap jet; |
| 35 | +% plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2), |
| 36 | +% legend(num2str(Zlist)), |
| 37 | +% set(gca, 'FontSize', 20); %for labels |
| 38 | +% xlabel('z'), |
| 39 | +% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr') |
| 40 | +% sum(MergerRateByRedshiftByZ(:,1)) %Total BBH merger rate at z=0 |
| 41 | +% [Zlist,MergerRateByRedshiftByZLanger]=CosmicHistoryIntegrator(filename,zlist,1); |
| 42 | +% [Zlist,MergerRateByRedshiftByZLambert]=CosmicHistoryIntegrator(filename,zlist,2); |
| 43 | +% plot(zlist, sum(MergerRateByRedshiftByZLambert,2), 'r', ... |
| 44 | +% zlist, sum(MergerRateByRedshiftByZLanger,2), 'b', 'LineWidth', 2), |
| 45 | +% set(gca, 'FontSize', 20), xlabel('z'), legend('Langer & Norman','Lambert') |
| 46 | +% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr'), |
| 47 | + |
| 48 | + |
| 49 | + |
| 50 | +%define constants |
| 51 | +global Mpcm; |
| 52 | +global Mpc; |
| 53 | +global yr; |
| 54 | +Mpcm=1*10^6 * 3.0856775807e16; %Mpc in meters |
| 55 | +c=299792458; %speed of light, m/s |
| 56 | +Mpc=Mpcm/c; %Gpc in seconds |
| 57 | +yr=3.15569e7; %year in seconds |
| 58 | + |
| 59 | +if (nargin<2) |
| 60 | + error('Not enough input arguments.'); |
| 61 | +end; |
| 62 | +if (nargin<4), makeplots=0; end; |
| 63 | + |
| 64 | +%cosmology calculator |
| 65 | +[zvec,tL]=Cosmology(); |
| 66 | +%load COMPAS data |
| 67 | +[M1,M2,Z,Tdelay]=DataRead(filename, max(tL)); |
| 68 | +%metallicity-specific SFR |
| 69 | +[Zlist,SFR,SFRfractionZ]=Metallicity(Z,zvec,metallicitychoice); |
| 70 | + |
| 71 | + |
| 72 | +%Consider the contribution of every simulated binary to the merger rate |
| 73 | +%in every redshift bin by considering when it would have to be formed to |
| 74 | +%merge at that redshift and normalizing by the relevant |
| 75 | +%metallicity-specific star formation rate |
| 76 | +dz=zvec(2)-zvec(1); |
| 77 | +tLmerge=tL(floor(zlist/(zvec(2)-zvec(1)))+1); |
| 78 | +MergerRateByRedshiftByZ=zeros(length(zlist),length(Zlist)); |
| 79 | +for(i=1:length(M1)), |
| 80 | + Zcounter=find(Zlist==Z(i)); |
| 81 | + zformindex=find(tL>=Tdelay(i),1); |
| 82 | + for(k=1:length(zlist)), |
| 83 | + zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1); |
| 84 | + MergerRateByRedshiftByZ(k,Zcounter)=... |
| 85 | + MergerRateByRedshiftByZ(k,Zcounter)+... |
| 86 | + SFR(zformindex)*SFRfractionZ(zformindex,Zcounter)/Msimulated; |
| 87 | + end; |
| 88 | +end; |
| 89 | + |
| 90 | +if(makeplots==1), %make a set of default plots |
| 91 | + MakePlots(M1,M2,Z,Tdelay,zvec,SFR,SFRfractionZ,... |
| 92 | + zlist,Zlist,MergerRateByRedshiftByZ); |
| 93 | +end; |
| 94 | + |
| 95 | +end %end of CosmicHistoryIntegrator |
| 96 | + |
| 97 | + |
| 98 | +%Load the data stored in COMPAS output format from a file |
| 99 | +%Select only binary black hole mergers of interest, and return the |
| 100 | +%component masses, metallicities, and star formation to merger delay times |
| 101 | +function [M1,M2,Z,Tdelay]=DataRead(filename, tHubble) |
| 102 | + global BH |
| 103 | + if(exist(filename, 'file')~=2), |
| 104 | + error('Input file does not exist'); |
| 105 | + end; |
| 106 | + data = importdata(filename, '\t', 2); |
| 107 | + |
| 108 | + colnames=data.textdata(2,:); |
| 109 | + Z1index=find(strcmp(colnames,'Metallicity1')); |
| 110 | + M1index=find(strcmp(colnames,'M1')); |
| 111 | + M2index=find(strcmp(colnames,'M2')); |
| 112 | + tcindex=find(strcmp(colnames,'tc')); |
| 113 | + tformindex=find(strcmp(colnames,'tform')); |
| 114 | + type1index=find(strcmp(colnames,'stellarType1')); |
| 115 | + type2index=find(strcmp(colnames,'stellarType2')); |
| 116 | + nonRLOF=data.data(:,find(strcmp(colnames,'RLOFSecondaryAfterCEE')))==0; |
| 117 | + Pessimistic=data.data(:,find(strcmp(colnames,'optimisticCEFlag')))==0; |
| 118 | + tc=data.data(:,tcindex); |
| 119 | + tform=data.data(:,tformindex); |
| 120 | + Ttotal=(tc+tform)*1e6; %convert Megayears to years |
| 121 | + %select only binaries where both members are black holes at the last |
| 122 | + %step, the total delay time is less than the age of the Universe, there |
| 123 | + %is no Roche-lobe overflow immediately after the common-envelope phase, |
| 124 | + %and the Pessimistic common-envelope conditions are met |
| 125 | + select=Ttotal<tHubble & nonRLOF & Pessimistic & ... |
| 126 | + data.data(:,type1index)==14 & data.data(:,type2index)==14; |
| 127 | + |
| 128 | + M1=data.data(select,M1index); |
| 129 | + M2=data.data(select,M2index); |
| 130 | + Z=data.data(select,Z1index); |
| 131 | + Tdelay=Ttotal(select); |
| 132 | +end %end of DataRead |
| 133 | + |
| 134 | +%Compute the star formation rate and lookback time (in years) |
| 135 | +%for an array of redshifts |
| 136 | +function [zvec,tL]=Cosmology() |
| 137 | + global Mpcm |
| 138 | + global Mpc |
| 139 | + global yr |
| 140 | + zmax=10; dz=0.001; zvec=0:dz:zmax; |
| 141 | + Nz=length(zvec); |
| 142 | + |
| 143 | + %Planck cosmology |
| 144 | + OmegaM=0.236+0.046; %2012arXiv1212.5226H |
| 145 | + OmegaL=0.718; |
| 146 | + Ho=100*0.697*1000/Mpcm; %in sec; H=69.7 |
| 147 | + Dh=1/Ho; |
| 148 | + E=sqrt(OmegaM.*(1+zvec).^3+OmegaL); %Hogg, astro-ph/9905116, Eq. 14 |
| 149 | + Dc=Dh*dz*cumsum(1./E); %Hogg, Eq. 15 |
| 150 | + Dm=Dc; %Hogg, Eq. 16, k=0; |
| 151 | + Dl=(1+zvec).*Dm; %Hogg, Eq. 20 |
| 152 | + %see also Eq. (1.5.46) in Weinberg, "Cosmology", 2008 |
| 153 | + dVc=4*pi*Dh^3*(OmegaM*(1+zvec).^3+OmegaL).^(-0.5).*(Dc/Dh).^2*dz/Mpc^3; |
| 154 | + Vc=cumsum(dVc); |
| 155 | + dtL=(1/Ho)*dz./(1+zvec)./E/yr; %lookback time, (30) of Hogg |
| 156 | + tL=cumsum(dtL); |
| 157 | +end %end of Cosmology |
| 158 | + |
| 159 | + |
| 160 | +%Compute the fraction of star formation that happens in a given metallicity |
| 161 | +%bin as a function of redshift |
| 162 | +function [Zlist,SFR,SFRfractionZ]=Metallicity(zvec,Z) |
| 163 | + %M_/odot per Mpc^3 per year -- Neijssel+ 2019 preferred model |
| 164 | + %would be SFR=0.015*(1+zvec).^2.7./(1+((1+zvec)/2.9).^5.6) in Madau & Dickinson, 2014, (15) |
| 165 | + SFR=0.01*(1+zvec).^2.77./(1+((1+zvec)/2.9).^4.7); |
| 166 | + Zmean=0.035.*10.^(-0.23*zvec); |
| 167 | + Zmu=log(Zmean)-0.39^2/2; |
| 168 | + dlogZ=0.01; |
| 169 | + logZvec=-12:dlogZ:0; %natural log |
| 170 | + dPdlogZ=1/0.39/sqrt(2*pi)*exp(-(logZvec'-Zmu).^2/2/0.39^2); |
| 171 | + dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise |
| 172 | + Zrange=log(max(Z))-log(min(Z)); %ugly correction for not including tails |
| 173 | + PdrawZ=1/Zrange; |
| 174 | + minlogZindex=find(exp(logZvec)>=min(Z),1, 'first'); |
| 175 | + maxlogZindex=find(exp(logZvec)<=max(Z),1, 'last'); |
| 176 | + dPdlogZ(minlogZindex,:)=dPdlogZ(minlogZindex,:)+sum(dPdlogZ(1:minlogZindex,:),1)*dlogZ/(sum(Z==min(Z))/length(Z))*PdrawZ; |
| 177 | + dPdlogZ(maxlogZindex,:)=dPdlogZ(maxlogZindex,:)+sum(dPdlogZ(maxlogZindex:end,:),1)*dlogZ/(sum(Z==max(Z))/length(Z))*PdrawZ; |
| 178 | + dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise |
| 179 | +end %end of Metallicity |
| 180 | + |
| 181 | + |
| 182 | +%Make a set of default plots |
| 183 | +function MakePlots(M1,M2,Z,Tdelay,zvec,SFR,SFRfractionZ,... |
| 184 | + zlist,Zlist,MergerRateByRedshiftByZ) |
| 185 | + |
| 186 | + figure(1),colormap jet; |
| 187 | + plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2), |
| 188 | + legend(num2str(Zlist)), |
| 189 | + set(gca, 'FontSize', 20); %for labels |
| 190 | + xlabel('z'), |
| 191 | + ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr') |
| 192 | + disp(['Total BBH merger rate at z=0: ', ... |
| 193 | + num2str(sum(MergerRateByRedshiftByZ(:,1))),... |
| 194 | + ' per Gpc^3 per year']); |
| 195 | + |
| 196 | + |
| 197 | + figure(2), colormap jet; |
| 198 | + scatter(M1,M2,20,log(Z)/log(10),'filled'); |
| 199 | + set(gca, 'FontSize', 20); %for labels |
| 200 | + H=colorbar; H.Label.String='log_{10} metallicity'; |
| 201 | + xlabel('Primary BH mass [M_o]'), ylabel('Secondary BH mass [M_o]'); |
| 202 | + |
| 203 | + figure(3), colormap jet; |
| 204 | + scatter(M1+M2,log(Tdelay/1e6)/log(10),20,log(Z)/log(10),'filled'); |
| 205 | + set(gca, 'FontSize', 20); %for labels |
| 206 | + H=colorbar; H.Label.String='log_{10} metallicity'; |
| 207 | + xlabel('Total BBH mass [M_o]'), ylabel('log(Tdelay/Myr)'); |
| 208 | + |
| 209 | + figure(4), colormap jet; |
| 210 | + plot(zvec, SFR) |
| 211 | + set(gca, 'FontSize', 20); %for labels |
| 212 | + xlabel('z'), ylabel('Star-formation rate, M_o per Gpc^3 per yr'); |
| 213 | + |
| 214 | + |
| 215 | + figure(5), colormap jet; |
| 216 | + plot(zvec, SFRfractionZ) |
| 217 | + set(gca, 'FontSize', 20); %for labels |
| 218 | + legend(num2str(Zlist)) |
| 219 | + xlabel('z'), ylabel('Z-specific SFR fraction'); |
| 220 | + |
| 221 | +end %end of MakePlots |
| 222 | + |
| 223 | + |
| 224 | + |
0 commit comments