Implementing Eddington limit for BHs based on spinup/ISCO (#570)

* changed BH eddington accretion to Bardeen 1970 prescription that adjusts the accretion based on the BH ISCO due to spinup since BH formation

* added one more BH eddington calculation for the newly added dme calculation

* fixed small typo in Marchant eddlim calc

* removed write statements I accidentally left in

* removed tb in a couple of the dme calls

* removed testing write statement

* removed dummy counter for marking when BH formed

Author: michaelzevin

cosmic/src/const_bse.h

@@ -14,8 +14,8 @@
      &               qcflag,eddlimflag,bhspinflag,aic,rejuvflag,
      &               htpmb,ST_cr,ST_tide,bdecayfac,grflag,
      &               bhms_coll_flag
-      REAL*8 don_lim,acc_lim
-      COMMON /MTVARS/ don_lim,acc_lim
+      REAL*8 don_lim,acc_lim,Mbh_initial
+      COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial
       INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn
       COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn
       INTEGER pisn_track(2)

cosmic/src/evolv2.f

@@ -185,6 +185,7 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
       REAL*8 bhspin(2)
       REAL*8 delet,delet1,dspint(2),djspint(2),djtx(2)
       REAL*8 dtj,djorb,djgr,djmb,djt,djtt,rmin,rdisk
+      REAL*8 etaBH,maxspinBH
 *
       INTEGER pulsar
       INTEGER mergemsp,merge_mem,notamerger,binstate,mergertype
@@ -270,6 +271,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
       ngtv2 = -2.d0
       twopi = 2.d0*ACOS(-1.d0)
 
+      Mbh_initial = 0.d0
+
 
 * disrupt tracks if system get disrupted by a SN during the common
 * envelope
@@ -523,6 +526,18 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
 * Just in case the wind mass loss rates are *very* high
 *
             dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(3-k)
+* For BHs, follow Marchant et al. 2017
+            if(kstar(3-k).eq.14)then
+               maxspinBH = 6.d0**(1.d0/2.d0) * Mbh_initial
+               if(mass(3-k).lt.maxspinBH)then
+                  etaBH = 1.d0 -
+     &    (1.d0 - (mass(3-k)/(3.d0*Mbh_initial))**(2.d0))**(1.d0/2.d0)
+               else
+                  etaBH = 0.42
+               endif
+               dme = 1.04e-3*eddfac*(1.d0/(1.d0 + zpars(11)))
+     &    *(1.d0/etaBH)*rad(3-k)
+            endif
 
 *
 * Calculate wind mass loss from the previous timestep.
@@ -1402,10 +1417,22 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
 *
 *
 * Force new NS or BH to have a birth spin peirod and magnetic field.
-* 
+*
          if(kstar(k).eq.13.or.kstar(k).eq.14)then
 * re-calculate the Eddington limit of the newly formed CO
-            dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(k)  
+            dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(k)
+* For BHs, follow Marchant et al. 2017
+            if(kstar(k).eq.14)then
+               maxspinBH = 6.d0**(1.d0/2.d0) * Mbh_initial
+               if(mass(k).lt.maxspinBH)then
+                  etaBH = 1.d0 -
+     &    (1.d0 - (mass(k)/(3.d0*Mbh_initial))**(2.d0))**(1.d0/2.d0)
+               else
+                  etaBH = 0.42
+               endif
+               dme = 1.04e-3*eddfac*(1.d0/(1.d0 + zpars(11)))
+     &    *(1.d0/etaBH)*rad(k)
+            endif
             if(tphys-epoch(k).lt.tiny)then
                if(kstar(k).eq.13.and.pulsar.gt.0)then
 *                  write(93,*)'birth start: ',tphys,k,B_0(k),ospin(k)
@@ -1906,7 +1933,19 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
 * Eddington limit for accretion on to the secondary in one orbit.
 *
  8    dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(j2)*tb
-      
+* For BHs, follow Marchant et al. 2017
+      if(kstar(j2).eq.14)then
+         maxspinBH = 6.d0**(1.d0/2.d0) * Mbh_initial
+         if(mass(j2).lt.maxspinBH)then
+            etaBH = 1.d0 -
+     &    (1.d0 - (mass(j2)/(3.d0*Mbh_initial))**(2.d0))**(1.d0/2.d0)
+         else
+            etaBH = 0.42
+         endif
+         dme = 1.04e-3*eddfac*(1.d0/(1.d0 + zpars(11)))
+     &    *(1.d0/etaBH)*rad(j2)*tb
+      endif
+
       supedd = .false.
       novae = .false.
       disk = .false.
@@ -2612,6 +2651,18 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
 * Just in case the wind mass loss rates are *very* high
 *
             dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(3-k)*tb
+* For BHs, follow Marchant et al. 2017
+            if(kstar(3-k).eq.14)then
+               maxspinBH = 6.d0**(1.d0/2.d0) * Mbh_initial
+               if(mass(3-k).lt.maxspinBH)then
+                  etaBH = 1.d0 -
+     &    (1.d0 - (mass(3-k)/(3.d0*Mbh_initial))**(2.d0))**(1.d0/2.d0)
+               else
+                  etaBH = 0.42
+               endif
+               dme = 1.04e-3*eddfac*(1.d0/(1.d0 + zpars(11)))
+     &    *(1.d0/etaBH)*rad(3-k)*tb
+            endif
 
             if(neta.gt.tiny)then
                if(beta.lt.0.d0)then !PK. following startrack

cosmic/src/hrdiag.f

@@ -800,6 +800,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
                         endif
                      endif
                      mt = mrem
+* Store the initial BH mass for calculating the ISCO later
+                     if(Mbh_initial.eq.0)then
+                        Mbh_initial = mt
+                     endif
                   endif
                endif
             endif
@@ -1176,7 +1180,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
                         endif
                      endif
                      mt = mrem
-
+* Store the initial BH mass for calculating the ISCO later
+                     if(Mbh_initial.eq.0)then
+                        Mbh_initial = mt
+                     endif
                      endif
                   endif
                endif
@@ -1271,6 +1278,10 @@ SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
 *        Black hole
 *
          mc = mt
+* Store the initial BH mass for calculating the ISCO later
+         if(Mbh_initial.eq.0)then
+            Mbh_initial = mt
+         endif
          lum = 1.0d-10
          r = 4.24d-06*mt
       endif