notebooks updates

notebooks/ShockPoweredTransients-Fang2021.py

@@ -54,6 +54,9 @@ import numpy as np
 # For a thin equatorial disk, $f_{\Omega} \sim h/r$.
 #
 
+# assumed f_omega in the following
+f_omega = 1.
+
 # ### Steady wind approximation
 #
 # Assume a steady wind with
@@ -66,7 +69,9 @@ import numpy as np
 # $
 #    n \simeq \dot M/(4\pi f_{\Omega}r^2 v_w m_p) = A / (m_p r^2)
 # $
+#
 # with
+#
 # $
 # A \equiv \dot M/(4\pi f_{\Omega} v_w) 
 # $
@@ -95,7 +100,7 @@ M_dot = np.logspace( -5., 0., 100)*u.M_sun/u.year
 plt.gcf().clear()
 plt.figure(figsize=(6,6))
 for z in v_w:
-    A = steady_wind_A(M_dot, z, 1.)
+    A = steady_wind_A(M_dot, z, f_omega)
     label_vw = ("v$_w=$%.0f km s$^{-1}$" % z.value)
     plt.plot(
         M_dot,
@@ -107,8 +112,10 @@ plt.yscale('log')
 plt.xscale('log')
 plt.legend()
 
-plt.title("steady wind loss parameter A")
+plt.title("steady wind loss parameter A ($f_{\Omega}=%.2f$)" % f_omega)
 plt.show()
+
+
 # -
 
 # ## Shocks
@@ -129,13 +136,140 @@ plt.show()
 # **Kinetic power of the forward shock**:
 #
 # $
-# L_{sh} = \frac{9\pi}{8} f_{\Omega} m_p n_{sh} v_{sh}^3 R_{sh}^2 = \frac{9}{32} \dot M \frac{v_{sh}^3}{v_w} = \frac{9\pi}{8} A f_{\Omega} v_{sh}^3
+# L_{sh} = \frac{9\pi}{8} f_{\Omega} m_p n_{sh} v_{sh}^3 R_{sh}^2 = \frac{9}{32} \dot M \frac{v_{sh}^3}{v_w} = \frac{9\pi}{8} A f_{\Omega} v_{sh}^3 \ \ (1)
+# $
+#
+# with $n_{sh}(R_{sh})$: upstream density ahead of shock
+#
+
+# +
+def kinetic_power(A, v_sh, fomega=1):
+    "kinetic power of forward shock"
+    return 9.*math.pi/8*A*fomega*v_sh**3
+
+
+v_sh = (np.logspace( 3, 5, 3))*u.km/u.s
+A = np.logspace( 11., 18., 100)*u.g/u.cm
+
+plt.gcf().clear()
+plt.figure(figsize=(8,8))
+for z in v_sh:
+    l_kin = kinetic_power(A, z, f_omega)
+    label_vw = ("v$_{sh}=$%.0f km s$^{-1}$" % z.value)
+    plt.plot(
+        A,
+        l_kin.to(u.erg/u.s),
+        label=label_vw)
+plt.xlabel("steady wind loss parameter A ($g/cm$)",fontsize=14)
+plt.ylabel("$L_{sh} (erg/s)$",fontsize=14)
+plt.yscale('log')
+plt.xscale('log')
+plt.legend()
+
+plt.title("Kinetic power of forward shock ($f_{\Omega}=%.2f$)" % f_omega)
+plt.show()
+# -
+
+# **Shock heating:**
+#
+# Gas immediatly behind the shock is heated to: $kT \simeq \frac{3}{16} \mu m_p v_{sh}^2 \approx 11 v_{8.5}^2$ keV.
+#
+# with $v_{sh}=v_{sh}/(3000 km s^{-1}$.
+
+# Mean molecular weight of the fully ionized gas (solar composition): $\mu=0.62$
+#
+# (would be $\mu=2$ if it would be hydrogen-poor gas)
+
+mu=0.62
+
+
+# +
+def shock_temperature(v_sh,mu):
+    "gas temperature behind the shock"
+    return 3./16.*mu*const.m_p*v_sh**2
+
+v_sh = (np.logspace( 3, 5, 100))*u.km/u.s
+T_sh = shock_temperature(v_sh, mu)
+
+plt.gcf().clear()
+plt.figure(figsize=(6,6))
+
+plt.plot(
+    v_sh,
+    T_sh.to(u.keV))
+
+plt.xlabel("shock speed ($km/s$)",fontsize=14)
+plt.ylabel("shock temperature ($keV$)",fontsize=14)
+plt.yscale('log')
+plt.xscale('log')
+
+plt.title("Temperature of gas behind shock ($\mu=%.2f$)" % mu)
+plt.show()
+
+# -
+
+# Large photoelectric opacity of the external medium: most of $L_{sh}$ is absorped and reprocessed via continuum and line emission into the optical wavelengths.
+#
+# Time delay of optical emission, as reprocessed emission from nearby the shock must propagate through the column of external medium: $\Sigma = \int_{R_sh}^{\infty} ndr \sim n_{sh} R_{sh}$.
+#
+# Optical diffussion time scale:
+#
+# $t_{diff} \approx \tau_{opt} \cdot (R_{sh}/c)$
+#
+# with the optical depth $\tau_{opt} \equiv \sum \sigma_{opt}$ ($\sigma_{opt}$ is the effective cross section at visual wavelengths).
+#
+# Adiabatic losses can be neglected as long as $t_{diff} \leq t_{dyn}$ with the expansion time scale of the shocked gas $t_{dyn} = R_{sh}/v_{sh}$:
+#
+# $
+#   \tau_{opt} \lesssim c/v_{sh} \ \ \ (3)
 # $
+
+# This condition is satisfied at the critical time
+#
+# $
+# t \gtrsim t_{pk} \approx \frac{c}{v_{sh}^2 n_{sh}} \sigma_{opt} = \frac{\dot M \kappa_{opt}}{4\pi f_{\Omega}c v_w} = \frac{A \kappa_{opt}}{c}  \ \ (4)
+# $
+#
+# with the optical opacity $\kappa_{opt}\equiv \sigma_{opt}/m_p$ (opacity: cross section for absorbing photons per units mass material; units [m$^2$/kg]).
 #
-# - $n_{sh}(R_{sh})$: upstream densitiy ahead of shock
+# Label $t_{pk}$ is introduced here, as the critical time often corresponds to the peak of the light curve.
 #
-# **Understand: derived from $1/2 m v^2$?**
+# Any diffusion through the eject is ignored in this step.
 
+# +
+def peak_time(a,kappa):
+    "peak time"
+    return a*kappa/const.c
+
+kappa_es = 0.38 *u.cm**2/u.g
+t_pk = peak_time(A, kappa_es)
+
+plt.gcf().clear()
+plt.figure(figsize=(6,6))
+
+plt.plot(
+    A,
+    t_pk.to(u.s))
+
+plt.xlabel("steady wind loss parameter A ($g/cm$)",fontsize=14)
+plt.ylabel("peak time $t_{pk}$ ($s$)",fontsize=14)
+plt.yscale('log')
+plt.xscale('log')
 
+plt.title("Peak time ($\kappa=%.2f$ cm$^2$/g )" % kappa_es.value)
+plt.show()
+# -
 
+# **Shock velocity is function of peak luminosity and peak time**:
+#
+# $
+# v_{sh} = \left( \frac{8}{9\pi} \frac{L_{pk} \kappa_{opt}}{ct_{pk}f_{\Omega}} \right)^{1/3}
+# $
+# (derived from combining (1) and (4))
+#
+# Assume here that $L_{pk} \approx L_{sh}$ (all optical emission is shock powered, without any contribution from e.g. radiactive decay).
 
+# ## Calorimetric Technique
+#
+#
+#