Should luminosity of Naked Helium stars be calculated using a fit for GB and AGB stars?? SSE paper does not really say.

[scottcoughlin2014]

TODO: add the example here

[astroJeff]

Nam found this minimal working example. When lambdaf is set to 0.5, the code breaks somewhere. When lambdaf=1, it runs smoothly.

It's hard to tell how serious of an issue this is, but when we get a chance, I suggest we work on debugging it.

from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.evolve import Evolve

m1 =6.242963173038605
m2 = 2.1269020737706663
porb =4900.260678899167
ecc = 0.8803069620359644
tphysf = 100
k1 = 1
k2 = 1
z = 0.02

single_binary = InitialBinaryTable.InitialBinaries(m1=m1,
                                                   m2=m2,
                                                   porb=porb,
                                                   ecc=ecc,
                                                   tphysf=tphysf,
                                                   kstar1=1,
                                                   kstar2=1,
                                                   metallicity=z)

BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.5, 'mxns': 2.5, 'beta': 0.125, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -1.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 2, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 0, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'zsun' : 0.02}

BSEDict1 = {
    'xi': 0.5,
    'bhflag': 1,
    'neta': 0.5,
    'windflag': 3,
    'wdflag': 0,
    'alpha1': 0.5,
    'pts1': 0.001,
    'pts3': 0.02,
    'pts2': 0.01,
    'epsnov': 0.001,
    'hewind': 1,
    'bwind': 0.0,
    'lambdaf': 0.5,
    'mxns': 2.5,
    'beta': -1,
    'tflag': 1,
    'acc2': 1.5,
    'nsflag': 3,
    'ceflag': 0,
    'eddfac': 1.0,
    'ifflag': 0,
    'sigma': 150.0,
    'gamma': -2.0,
    'pisn': 45.0,
    'bhsigmafrac': 1.0,
    'polar_kick_angle': 90,
    'qcrit_array': [
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        0.0, 0.0
    ],
    'cekickflag': 2,
    'cehestarflag': 0,
    'cemergeflag': 0,
    'ecsn': 2.5,
    'ecsn_mlow': 1.4,
    'aic': 1,
    'ussn': 0,
    'sigmadiv': 9999999999999,
    'qcflag': 2,
    'eddlimflag': 0,
    'fprimc_array': [
        2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0,
        2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0,
        2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0, 2.0 / 21.0
    ],
    'bhspinflag': 0,
    'bhspinmag': 0.0,
    'rejuv_fac': 1.0,
    'rejuvflag': 0,
    'htpmb': 1,
    'ST_cr': 1,
    'ST_tide': 1,
    'bdecayfac': 1,
}

for k,v in BSEDict1.items():
    BSEDict[k] = v

bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)

print(bpp)

[scottcoughlin2014]

Eventually we should change this issue because it actually is not rooted in lambaf value at all but here we go.

The issue here is that

Below is the line that returns a NaN result from hrdiag.f

 816 *
 817 * Helium Shell Burning
 818 *
 819             kw = 8
 820             lum = lgbtf(aj,GB(8),GB,tscls(4),tscls(5),tscls(6))
 821             r = rhehgf(mt,lum,rx,lums(2))
 822             rg = rhegbf(lum)
 823             if(r.ge.rg)then
 824                kw = 9
 825                r = rg
 826             endif
 827             mc = mcgbf(lum,GB,lums(6))
 828             mtc = MIN(mt,1.45d0*mt-0.31d0)
 829             mcmax = MIN(mtc,MAX(mch,0.773d0*mass-0.35d0))

lgbtf returns a NaN result because the output of tinf2-t is negative and we take a root of a negative which is bad. This implies the time scales used for Naked Helium Stars (which in the SSE paper is brushed under the rug and they simple follow the logic of the GB (for those equations continue reading). The function lgbtf is from zfuncs.f and is as follows

 583       real*8 FUNCTION lgbtf(t,A,GB,tinf1,tinf2,tx)
 584       implicit none
 585       real*8 t,A,GB(10),tinf1,tinf2,tx
 586 *
 587 * A function to evaluate L given t for GB, AGB and NHe stars
 588 *
 589       if(t.le.tx)then
 590          lgbtf = GB(4)*(((GB(5)-1.d0)*A*GB(4)*(tinf1 - t))**
 591      &                              (GB(5)/(1.d0-GB(5))))
 592       else
 593          lgbtf = GB(3)*(((GB(6)-1.d0)*A*GB(3)*(tinf2 - t))**
 594      &                              (GB(6)/(1.d0-GB(6))))
 595       endif
 596 *
 597       return
 598       end

tx=tscls(6) tinf1=tscls(4) tinf2=tscls(5) are all set in the code below (see screen shot of SSE paper for the relevant equations in the paper. (this is from star.f)

114 * Set GB timescales
115          tscls(4) = tscls(1) + (1.d0/((GB(5)-1.d0)*GB(1)*GB(4)))*
116      &              ((GB(4)/lums(3))**((GB(5)-1.d0)/GB(5)))
117          tscls(6) = tscls(4) - (tscls(4) - tscls(1))*((lums(3)/lums(6))
118      &              **((GB(5)-1.d0)/GB(5)))
119          tscls(5) = tscls(6) + (1.d0/((GB(6)-1.d0)*GB(1)*GB(3)))*
120      &              ((GB(3)/lums(6))**((GB(6)-1.d0)/GB(6)))