diff --git a/examples/Params.ini b/examples/Params.ini index 0e9282dba..3b5dc6f56 100644 --- a/examples/Params.ini +++ b/examples/Params.ini @@ -419,6 +419,15 @@ alpha1 = [1.0, 1.0] ; Default: 0.0 lambdaf = 0.0 +; ce2stageflag +; Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022 +; +; Options: +; 0 - Use the standard energy formalism +; 1 - Use the two-stage CE formalism from Hirai & Mandel 2022 +; Default: 0 +ce2stageflag = 0 + ; ceflag ; Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses) ; diff --git a/src/cosmic/consts.py b/src/cosmic/consts.py index f3d94c627..d2967c9c8 100644 --- a/src/cosmic/consts.py +++ b/src/cosmic/consts.py @@ -42,7 +42,7 @@ "epsnov", "eddfac", "gamma", "LBV_flag", ], "cevars": ["alpha1", "lambdaf", "qcrit_array"], - "ceflags": ["ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"], + "ceflags": ["ce2stageflag", "ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"], "flags": [ "tflag", "ifflag", "wdflag", "rtmsflag", "bhflag", "remnantflag", "maltsev_mode", "grflag", "bhms_coll_flag", "wd_mass_lim", "aic", @@ -63,4 +63,4 @@ "mixvars": ["rejuv_fac"], "metvars": ["zsun"], "rand1": ["randomseed"] -} \ No newline at end of file +} diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index fca181857..59ada51f9 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -865,6 +865,23 @@ } ] }, + { + "name": "ce2stageflag", + "description": "Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022", + "type": "dropdown", + "options-preface": "", + "options": [ + { + "name": 0, + "description": "Use the standard energy formalism", + "default": true + }, + { + "name": 1, + "description": "Use the two-stage CE formalism from Hirai & Mandel 2022" + } + ] + }, { "name": "ceflag", "description": "Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses)", diff --git a/src/cosmic/src/SSE/SSE_zfuncs.f b/src/cosmic/src/SSE/SSE_zfuncs.f index 25b78578a..7bce39384 100644 --- a/src/cosmic/src/SSE/SSE_zfuncs.f +++ b/src/cosmic/src/SSE/SSE_zfuncs.f @@ -845,6 +845,94 @@ real*8 FUNCTION lmcgbf(mc,GB) * return end +*** + real*8 FUNCTION tonset(tmin,z) + implicit none + real*8 tmin,z +* +* A function to evaluate Tonset (Teff when Mconv reaches 1% of its max) +* from Mandel, Hirai & Picker 2024 +* + tonset = tmin /min(0.695d0-0.057d0*log10(z),0.95d0) +* + return + end +*** + real*8 FUNCTION menvmax(kw,m,z) + implicit none + integer kw + real*8 m,z + real*8 mcagbf,mifmcf + external mcagbf +* +* A function to evaluate the maximum mass of the convective envelope +* from Picker, Hirai & Mandel 2024 +* +* I'm not sure how HG stars should be treated ... +* + if(kw.ge.10)then + menvmax = 0.0d0 + elseif(kw.lt.2)then + if(m.gt.1.25d0)then + menvmax = 0.0d0 + else + menvmax = m + if(m.gt.0.35d0)then + menvmax = 0.35d0*(1.25d0-m)*(1.25d0-m)/0.81d0 + endif + endif + else + mifmcf = -0.023*log10(z)-0.0023 + menvmax = max(m-mcagbf(m)*(1.0d0+mifmcf), 0.0d0) + endif +* + return + end +*** + real*8 FUNCTION mconvenv(kw,m,z,teff,tmin,tset,aj,tm) + implicit none + integer kw + real*8 m,z,teff,tmin,tset,aj,tm + real*8 menvmax + external menvmax +* +* A function to evaluate the mass of the convective envelope +* from Picker, Hirai & Mandel 2024 +* + if(kw.ge.10)then + mconvenv = 0.0d0 + elseif(kw.lt.2)then + if(m.gt.1.25d0)then + mconvenv = 0.0d0 + else + mconvenv = menvmax(kw,m,z)*sqrt(sqrt(1.0d0-aj/tm)) + endif + else + mconvenv = menvmax(kw,m,z)/(1.0d0+exp(4.6d0* + & (tmin+tset-2.0d0*teff)/(tmin-tset))) + endif +* + return + end +*** + real*8 FUNCTION celamhe(m,z,mconv,mconvmax) + implicit none + real*8 m,z,mconv,mconvmax + real*8 m2,b1 +* +* A function to evaluate the binding energy factor of the convective envelope +* from Picker, Hirai & Mandel 2024 +* + m2 = 0.0023d0*log10(z)**2+0.0088d0*log10(z)+0.013d0 + b1 = m2*m-0.23d0 + if((mconv/mconvmax).gt.0.3d0)then + celamhe = exp(0.42*(mconv/mconvmax)+b1) + else + celamhe = exp(0.3*0.42+b1) + endif +* + return + end *** real*8 FUNCTION lHeIf(m,mhefl) implicit none diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 849b0fe06..748efc802 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -75,6 +75,13 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, REAL*8 rad(2),tms(2),lumin(2),B_0(2),bacc(2),tacc(2),epoch(2) REAL*8 menv_bpp(2),renv_bpp(2) REAL*8 ALPHA_CE + + REAL*8 met,reagb,tmin,tset,mconvmax + REAL*8 ragbf,menvmax,mconvenv,CELAMHE,tonset + EXTERNAL ragbf,menvmax,mconvenv,CELAMHE,tonset + REAL*8 mcgbtf,lmcgbf + EXTERNAL mcgbtf,lmcgbf + REAL*8 m1endstage1,m2endstage1 * * Initialize * @@ -120,6 +127,44 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RZAMS = RZAMSF(M01) LAMB1 = CELAMF(KW,M01,L1,R1,RZAMS,MENVD,LAMBDAF) ENDIF +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + IF (using_SSE.eq.1) THEN +* Estimate the mass of the convective envelope +* from the fits in Picker, Hirai, & Mandel 2024 + met = 10**(LOG10(zpars(14))/0.4) + reagb = ragbf(M1,LUMS(7),zpars(2)) + tmin = 1000.d0*((1130.d0*LUMS(7)/ + & (reagb**2.d0))**(1.d0/4.d0)) + tset = tonset(tmin,met) + teff1 = 1000.d0*((1130.d0*lumin(1)/ + & (rad1_bpp**2.d0))**(1.d0/4.d0)) + MENV = mconvenv(KW,M1,met,teff1,tmin,tset,AJ1,TM1) + mconvmax = menvmax(KW,M1,met) + LAMB1 = CELAMHE(M1,met,MENV,mconvmax) + ENDIF +* if > 8 Msun, the envelope participating in the first stage is only the convective one +* if < 2 Msun, the entire envelope participates in the first stage +* linear interpolation in between + IF(M1.GE.8.0d0)THEN + m1endstage1 = M1-MENV + ELSEIF(M1.LT.2.0d0)THEN + m1endstage1 = MC1 + ELSE + m1endstage1 = MC1 + (M1-MC1 - MENV) * (M1-2.0d0) /6.0d0 + ENDIF +* +* Calculate the binding energy of the convective envelope (multiplied by lambda). +* + EBINDI = M1*(M1-m1endstage1)/(LAMB1*R1) + ELSE +* +* Calculate the binding energy of the giant envelope (multiplied by lambda). +* + EBINDI = M1*(M1-MC1)/(LAMB1*R1) + ENDIF KW = KW2 CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,dtm,star2) CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS, @@ -127,13 +172,10 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & bhspin2,star2) OSPIN2 = JSPIN2/(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2) * -* Calculate the binding energy of the giant envelope (multiplied by lambda). -* - EBINDI = M1*(M1-MC1)/(LAMB1*R1) -* * If the secondary star is also giant-like add its envelopes energy. * Determine EORBI based on CEFLAG (CEFLAG=1 for de Kool prescription) * + m2endstage1 = M2 IF(KW2.GE.2.AND.KW2.LE.9.AND.KW2.NE.7)THEN MENVD = MENV/(M2-MC2) IF (using_METISSE.eq.1) THEN @@ -142,7 +184,44 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RZAMS = RZAMSF(M02) LAMB2 = CELAMF(KW,M02,L2,R2,RZAMS,MENVD,LAMBDAF) ENDIF - EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2) +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + IF (using_SSE.eq.1) THEN +* Estimate the mass of the convective envelope +* from the fits in Picker, Hirai, & Mandel 2024 + met = 10**(LOG10(zpars(14))/0.4) + reagb = ragbf(M2,LUMS(7),zpars(2)) + tmin = 1000.d0*((1130.d0*LUMS(7)/ + & (reagb**2.d0))**(1.d0/4.d0)) + tset = tonset(tmin,met) + teff2 = 1000.d0*((1130.d0*lumin(2)/ + & (rad2_bpp**2.d0))**(1.d0/4.d0)) + MENV = mconvenv(KW,M2,met,teff2,tmin,tset,AJ2,TM2) + mconvmax = menvmax(KW,M2,met) + LAMB2 = CELAMHE(M2,met,MENV,mconvmax) + ENDIF +* if > 8 Msun, the envelope participating in the first stage is only the convective one +* if < 2 Msun, the entire envelope participates in the first stage +* linear interpolation in between + IF(M2.GE.8.0d0)THEN + m2endstage1 = M2-MENV + ELSEIF(M1.LT.2.0d0)THEN + m2endstage1 = MC2 + ELSE + m2endstage1 = MC2 + (M2-MC2 - MENV) * (M2-2.0d0) /6.0d0 + ENDIF +* +* Calculate the binding energy of the convective envelope (multiplied by lambda). +* + EBINDI = EBINDI + M2*(M2-m2endstage1)/(LAMB2*R2) + ELSE +* +* Calculate the binding energy of the giant envelope (multiplied by lambda). +* + EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2) + ENDIF * * Calculate the initial orbital energy * @@ -167,7 +246,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, ENDIF EORBF = EORBI + EBINDI/ALPHA_CE * -* If the secondary is on the main sequence see if it fills its Roche lobe. +* If the secondary is on the main sequence. * IF(KW2.LE.1.OR.KW2.EQ.7)THEN SEPF = MC1*M2/(2.D0*EORBF) @@ -176,6 +255,63 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RL1 = RL(Q1) RL2 = RL(Q2) * +* If degenerate or giant secondary. +* + ELSE + SEPF = MC1*MC2/(2.D0*EORBF) + Q1 = MC1/MC2 + Q2 = 1.D0/Q1 + RL1 = RL(Q1) + RL2 = RL(Q2) + ENDIF +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + SEPF = m1endstage1*m2endstage1/(2.D0*EORBF) +* log end of the first stage with evolve_type=3 + q1_bpp = m1endstage1/m2endstage1 + q2_bpp = 1.d0/q1_bpp + rrl1_bpp = RC1/(RL(q1_bpp)*SEP_postCE) + rrl2_bpp = R2/(RL(q2_bpp)*SEP_postCE) + TB = (SEPF/AURSUN)* + & SQRT(SEPF/(AURSUN*(m1endstage1+m2endstage1))) + CALL writetab(jp,tphys,3.d0, + & m1endstage1,m2endstage1,KW1,KW2, + & SEPF,TB,ECC, + & rrl1_bpp,rrl2_bpp, + & AJ1,AJ2,tms1_bpp,tms2_bpp, + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M02,M01,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2), + & renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') +* +* Second stage: stable mass transfer of the radiative intershell +* Binary hardening formula from Picker, Hirai & Mandel 2024 +* + IF(m1endstage1.GT.MC1)THEN + SEPF = SEPF*((m1endstage1+m2endstage1)/(MC1+m2endstage1)) + & *(m1endstage1/MC1)**2 + & *EXP(-2*(m1endstage1-MC1)/m2endstage1) + EORBF = MC1*m2endstage1/(2.D0*SEPF) + ENDIF + Q1 = MC1/m2endstage1 + Q2 = 1.D0/Q1 + RL1 = RL(Q1) + RL2 = RL(Q2) + ENDIF +* +* If the secondary is on the main sequence see if it fills its Roche lobe. +* + IF(KW2.LE.1.OR.KW2.EQ.7)THEN +* * If cemergeflag is set, cause kstars without clear core-envelope * structure to merge automatically if they enter a CE * diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index 3b3314f9d..022a8faed 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -18,7 +18,9 @@ REAL*8 don_lim,acc_lim(2),Mbh_initial,smt_periastron_check COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial,smt_periastron_check + INTEGER ce2stageflag INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn + COMMON /CEFLAGS/ ce2stageflag COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn INTEGER pisn_track(2) COMMON /TRACKERS/ pisn_track diff --git a/src/cosmic/src/test_bse.f b/src/cosmic/src/test_bse.f index 983d70cb0..47e2a6561 100644 --- a/src/cosmic/src/test_bse.f +++ b/src/cosmic/src/test_bse.f @@ -123,6 +123,7 @@ PROGRAM bse WRITE(*,*)mass0(1),mass0(2),tphysf,tb,kstar(1),kstar(2),z,ecc WRITE(*,*)neta,bwind,hewind,alpha1(1),alpha1(2),lambdaf,windflag + WRITE(*,*)ce2stageflag WRITE(*,*)rtmsflag,ceflag,tflag,ifflag,wdflag,bhflag,remnantflag WRITE(*,*)mxns,idum,pts1,pts2,pts3 WRITE(*,*)sigma,beta,xi,acc2,epsnov,eddfac,gamma diff --git a/src/cosmic/tests/data/Params.ini b/src/cosmic/tests/data/Params.ini index 9aba74431..8c0078ce0 100644 --- a/src/cosmic/tests/data/Params.ini +++ b/src/cosmic/tests/data/Params.ini @@ -194,6 +194,11 @@ alpha1 = [1.0, 1.0] ; default=0.5 lambdaf = 0.5 +; ce2stageflag=0 uses the standard energy formalism for CEE +; ce2stageflag=1 uses the two-stage formalism CEE from Hirai & Mandel 2022 +; default=0 +ce2stageflag=1 + ; ceflag=1 used the method from de Kool 1990 for setting the initial orbital energy ; ceflag=0 does not use this method (uses the core mass to calculate initial orbital energy) ; default=0 diff --git a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 index eb582d7ab..9a5e289ff 100644 Binary files a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 and b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 differ diff --git a/src/cosmic/tests/data/kick_initial_conditions.h5 b/src/cosmic/tests/data/kick_initial_conditions.h5 index 0323919eb..10aeee830 100644 Binary files a/src/cosmic/tests/data/kick_initial_conditions.h5 and b/src/cosmic/tests/data/kick_initial_conditions.h5 differ