-->

Wednesday, November 21, 2012

test


Using the numbers we measure for $L(H\alpha)$, I'll re-derive key quantities about SN 2009ip.


This analysis implicitly assumes that H$\alpha$ is optically thick at $10^4$ K. I believe we should make that assumption explicit in the text.





In [1]:


from astropy import units, constants
pc = units.pc.in_units('cm')
Iha = 2.9e-14 # erg s^-1 cm^-2 A^-1
D = 24e6 * pc
Lha = Iha * D**2
print Iha," erg/s/cm^2/A"
print D," cm"
print Lha," erg/s/A"









2.9e-14  erg/s/cm^2/A
7.40562619552e+25 cm
1.59045568109e+38 erg/s/A








In [2]:


# Use the wavelength form of the Planck function
h = constants.cgs.h
k = constants.cgs.k_B
c = constants.cgs.c
def planck(lam, temp):
""" erg s^-1 cm^-2 A^-1 sr^-1 """
return 2*h*c**2 / lam**5 * (exp(h*c/(k*temperature*lam)) - 1)**-1







In [3]:


angstrom = units.angstrom.in_units('cm') # angstroms
wave = 6562.8*angstrom
temperature = 1e4 # K
print "sanity check: wave * 1e8 = 6562.8 = ",wave*1e8
print "planck(6562.8,1e4) = ",planck(wave,temperature)," erg/s/cm^2/cm"
print "planck(6562.8,1e4) = ",planck(wave,temperature)/1e8," erg/s/cm^2/A"









sanity check: wave * 1e8 = 6562.8 =  6562.8
planck(6562.8,1e4) = 1.22966749446e+15 erg/s/cm^2/cm
planck(6562.8,1e4) = 12296674.9446 erg/s/cm^2/A







We're trying to solve an equation of the form:
$$ L(H\alpha) = \pi R^2 B_\lambda(T) d\lambda $$
Since our luminosity is in erg s$^{-1}$ $\dot{A}^{-1}$, we don't have to worry about the bandwidth.

$$ R = (L(H\alpha) / \pi / B_\lambda(T))^{1/2}$$





In [4]:


Aemitting = Lha / (planck(wave,temperature)*pi/1e8)
print "Emitting area = ",Aemitting," cm^2"
R = (Aemitting/pi)**0.5
print "R = %e" % R, " cm"
au = units.AU.in_units('cm')
print "R = ",R/au," AU"
print "R = ",R/pc, " pc"









Emitting area =  4.11702975893e+30  cm^2
R = 1.144767e+15 cm
R = 76.5229413217 AU
R = 0.000370993688718 pc








In [5]:


# Light crossing time
tc = R / constants.cgs.c
day = units.day.in_units('s')
print "Light Crossing time at %i AU = %f days" % (R/au, tc/day)









Light Crossing time at 76 AU = 0.441960 days








In [6]:


# What if we assume it's a sphere instead of a disk?
# area -> 4 * pi * r^2
R = (Aemitting/(4*pi))**0.5
print "R = %e" % R, " cm"
au = units.AU.in_units('cm')
print "R = ",R/au," AU"
print "R = ",R/pc, " pc"
tc = R / constants.cgs.c
print "Light Crossing time at %i AU = %f days" % (R/au, tc/day)









R = 5.723835e+14  cm
R = 38.2614706609 AU
R = 0.000185496844359 pc
Light Crossing time at 38 AU = 0.220980 days







Because the brightening was sudden, and because a uniform disk is somewhat implausible, we should come up with a different measurement for a disk with an inner hole. I'll use fiducial numbers of 13000 $km s^{-1}$ and $\Delta t \sim 40$ days.





In [14]:


kms = (units.km/units.s).in_units('cm/s')
vejecta = 13e3 * kms
deltat = 40 * day
r_collision = vejecta*deltat
print "Inner disk radius = ",r_collision," or ",r_collision/au," AU"
print "Light crossing time at inner disk radius = ",r_collision / constants.cgs.c / day, "days"
H = Aemitting/r_collision/2/pi
print "Height needed if this is ALL coming from the inner part of a disk: ",H," cm or ",H/au," AU"
print "The speed required to achieve that height in 2 days is ", H / (2*day) / 1e5, "km/s"









Inner disk radius =  4.4928e+15  or  300.325130229  AU
Light crossing time at inner disk radius = 1.73453329503 days
Height needed if this is ALL coming from the inner part of a disk: 1.45843491139e+14 cm or 9.74903522735 AU
The speed required to achieve that height in 2 days is 8440.01684832 km/s







We then want $$\int_{r_{inner}}^{r_{outer}} 2 \pi r dr = A_{emitting}$$ or $$ r_{outer}^2 - r_{inner}^2 = A_{emitting} / \pi$$





In [8]:


router = (Aemitting/pi+r_collision**2)**0.5
print "Outer radius = ",router," or ",router/au," AU"
print "Delta-R = ",router-r_collision," or ",(router-r_collision)/au," AU"
tcross = (router-r_collision) / vejecta
print "The shock crossing time from inner->outer at %g km/s = %g days" % (vejecta/kms, tcross/day)









Outer radius =  4.6363501932e+15  or  309.920867958  AU
Delta-R = 1.43550193199e+14 or 9.59573772859 AU
The shock crossing time from inner->outer at 13000 km/s = 1.27805 days







I think this effectively leaves us with the conclusion that there is a 30 AU radius torus at 300 AU from the central source. That strikes me as somewhat implausible; 30 AU is enormous, particularly at densities $\sim10^{13}$. I'll try computing the mass in the disk assuming it's quite thin, $\sim0.1AU$.





In [9]:


height = 0.1 * au
volume = Aemitting*height
density = 1e13
nparticles = volume * density
msun = units.msolMass.in_units('g')
mass = nparticles * constants.cgs.m_p * 1.4
print "mass = ",mass,"g = ",mass/msun,"msun"









mass =  1.44223224375e+32 g =  72.5067741062 msun







That's ridiculous. I'll invert the assumption: say the disk/torus is 1 msun; how thin must it be?





In [10]:


mass = 1 * msun
nparticles = mass / (constants.cgs.m_p*1.4)
volume = nparticles/density
height = volume/Aemitting
print "height of disk = ",height,"cm = ",height/au," AU"









height of disk =  20632261267.2 cm =  0.0013791814797  AU







These numbers don't feel right to me. Why such an ultrathin disk? It would HAVE to puff up; the pressures are enormous.


What about the hypothesis that the brightening results from optically thin emission? I'll work under the assumption that H$\alpha$ is "just barely" optically thick (i.e. I'll cheat and use the optically thin approach);
in order to get the high H$\beta$ / H$\alpha$ ratio we actually need it to be pretty optically thick and under the Case B regime (Case C just makes it harder to get the observed $\alpha$/$\beta$ ratio).





In [18]:


density = 1e6 # assume pretty dense, but still plausibly case B
alphaB = 2.6e-13 # cm^3 s^-1
Eha = h * c / wave
nphot_perang = Lha / Eha # phot/s/A
nphot = nphot_perang * 1 # phot/s
# solve for path length
# nphot = density**2 * alphaB * area * path_length
path_length = nphot / (density**2*alphaB*Aemitting)
print "Emitting path length for density = %0.1e" % density,"cm^-3 is ",path_length," = ",path_length/pc," pc"
# use the rise time of ~2 days as an upper limit on the emitting size scale
maxsize = 2*day*c
print "This is %i times larger than the largest allowed scale" % (path_length/maxsize)
print "Densities must be at least %e for optically thin emission to be plausible" % (density*sqrt(path_length/maxsize))









Emitting path length for density = 1.0e+06 cm^-3 is  4.90881531249e+19  =  15.9083870006  pc
This is 9475 times larger than the largest allowed scale
Densities must be at least 9.734331e+07 for optically thin emission to be plausible







So at densities $>5\times10^7$, if H$\alpha$ starts becoming optically thick but H$\beta$ remains optically thin, we could explain the brightening by a filled sphere of just barely optically thick emission.




Despite all this, we still don't have a great explanation for why the P Cygni profile disappears. We should revisit that if possible.





In [22]:


# limb-brightened shell model
half_sph_area = 2*pi*r_collision**2
print "Fraction of sphere that must be emitting: ",Aemitting/half_sph_area









Fraction of sphere that must be emitting:  0.0324616032628









In [1]:


from IPython.core.display import HTML

def blackbg():
"""Black Background code mirror theme."""
html = """
<style type="text/css">
.cm-s-ipython { background-color: black; color: lightblue; }
.cm-s-ipython span.cm-keyword {color: #00ff00; font-weight: bold;}
.cm-s-ipython span.cm-number {color: #ee88ee;}
.cm-s-ipython span.cm-operator {color: lime; font-weight: bold;}
.cm-s-ipython span.cm-meta {color: white;}
.cm-s-ipython span.cm-comment {color: cyan; font-style: italic;}
.cm-s-ipython span.cm-string {color: red;}
.cm-s-ipython span.cm-error {color: darkred;}
.cm-s-ipython span.cm-builtin {color: pink; font-weight: bold;}
.CodeMirror pre.CodeMirror-cursor {color: white; border-left: 1px solid white;}
.cm-s-ipython span.cm-variable {color: white;}
div.text_cell_input { background-color: black; color: white;}
div.text_cell { background-color: black; color: white;}
div#notebook { background-color: black; background-image: none; color: white;}
div.code_cell { background-color: black; color: white;}
div.metaedit .maintoolbar{ background-color: black; color: white;}
body{ background-color: black; color: white;}
#notebook_name { background-color: #333333; color: white; background-image:none;}
span.ui-widget-content{ background-color: #333333; color: white; background-image:none;}
div.ui-widget-content{ background-color: #333333; color: white; background-image:none;}
div#ui-button ui-widget ui-state-default ui-corner-all ui-button-text-only output_collapsed vbox{
background-color: black; color: white;}
div.input_area{ background-color: black; color: white;}
div.ui-widget-content { background-color: #333333; background-image: none;}
div.input_prompt {color:#6666CC; font-weight: bold;}
div.output_prompt {color:#CC6666; font-weight: bold;}
div.output_text {color:white;}
div.text_cell_render {color:white;}
span.item_name {color:white;}
.ui-widget-content a {color:white;}
div.tooltip .ui-corner-all .tooltiptext .smalltooltip {background-color: black; color:white;}
.pln {color:white}
.typ {color:#f9f}
.lit {color:#BBB}
.kwd {color:#3e3}
.metaedit {background: transparent; border-color: transparent; color:white;}
.ui-button-text {background: #112; color:white; border-color:666;}
.ui-menubar {background-color: transparent; background: transparent; background-image:none;}
.ui-widget-header {background-color: transparent; background: transparent; background-image:none;}
.ui-helper-clearfix {background-color: transparent; background: transparent; background-image:none;}
.ui-corner-all {background: #112; color: white;}
.ui-widget {background: #112; color: white;}
.ui-widget-content {background: #112; color: white;}
div.highlight {background:black;} /* needed for nbviewer */
</style>"""
return HTML(html)

blackbg()







Out[1]:









No comments: