Commit 68f06b06 authored by Gernot Maier's avatar Gernot Maier
Browse files

cleanup

parent a445141e
# ---
# jupyter:
# jupytext:
# formats: ipynb,md,py:light
# text_representation:
# extension: .py
# format_name: light
# format_version: '1.5'
# jupytext_version: 1.13.0
# kernelspec:
# display_name: Python 3 (ipykernel)
# language: python
# name: python3
# ---
# # Shock powered transients
#
# - [Fang, K. et al (2021)](https://arxiv.org/abs/2007.15742)
#
# Step throught this paper to capture the relevant points.
import astropy.units as u
import astropy.constants as const
import math
import matplotlib.pyplot as plt
import numpy as np
# #### Fang et al (2021) - Figure 1
#
# <img src="./figures/Fang2021-Figure1.png" alt=drawing width=500/>
# ## Section 2.1: Shock Dynamics and Thermal Expansion
#
# - spherical expanding homologoues ejecta of average velicity $\bar v_{ej}$
# - inner layers slower than other layers: $v_{ej}\propto r$ (homologoues velocity profile)
#
# ### External medium
#
# Density of external medium
#
# $
# n \equiv \rho / m_p
# $
# ($\rho$ = mass densitiy)
#
# with a radial profile of
#
# $
# n \propto r^{-k}
# $
# with $k\geq 2$
#
# The external medium is concentrated in a fractional solid angle $f_{\Omega}\leq 1$.
# 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
#
# - mass-loss rate $\dot M$
# - wind velocity $v_w$
#
# the density becomes
#
# $
# 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)
# $
#
# Assume here $k=2$, otherwise A(r) is a function of radius.
# #### Typical values
#
# From modeling interacting supernovae Smith (2014):
# - $\dot M \sim 10^{-4} - 1$ M$_{\odot}$ yr$^{-1}$
# - $v_w \sim 100 - 1000$ km s$^{-1}$
#
# (note the typo units in the Fang et al paper: A$_*$ should be in g/cm and not g/cm$^2$)
# +
# canonical value
A_star = 5.e11 *u.g/u.cm**2
def steady_wind_A(mdot, vw, fomega=1):
"mass loss A"
return mdot/(4.*math.pi*fomega*vw)
v_w = (np.linspace( 100, 1000, 5))*u.km/u.s
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, f_omega)
label_vw = ("v$_w=$%.0f km s$^{-1}$" % z.value)
plt.plot(
M_dot,
A.to(u.g/u.cm),
label=label_vw)
plt.xlabel("mass loss rate $\dot M$ (M$_{\odot}$/yr)",fontsize=14)
plt.ylabel("A (g/cm)",fontsize=14)
plt.yscale('log')
plt.xscale('log')
plt.legend()
plt.title("steady wind loss parameter A ($f_{\Omega}=%.2f$)" % f_omega)
plt.show()
# -
# ## Shocks
#
# Collision drives two shocks
# - forward shock into external medium
# - reverse shock into ejecta
#
# Shocks are radiative with rapid cooling of the gas behind both shocks --> gas accumulate in a thin shell which propagates outwards into the external medium.
#
#
# Shocks reach radius $R_{sh}\approx v_{sh}t$ by the time $t$ after the explosion.
# Shock speed $v_{sh}$ will reach ejecta speed $v_{ej}$ ([Metzer & Pejach (2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.471.3200M/abstract) **READ**), which reduces the power of the reverse shock with relation to the forward shock (**???**)
#
# $\rightarrow$ emission dominated by forward shock
# **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 \ \ (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 velocity ($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 times later than $t$:
#
# $
# 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]).
#
# Label $t_{pk}$ is introduced here, as the critical time often corresponds to the peak of the light curve.
#
# Any diffusion through the ejecta and corrections to $t_{diff}$ due its non-spherical geometry are 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).
# ## Calorimetry
#
# Conditions on the optical depth are very similar to the required conditions for particle acceleration and that the shock discontinuity is mediated by collisionless plasma. At times earlier than $t_{pk}$, the trapped radiations thickens the shock transition to macroscopic scale and does not allow any particle injection.
#
# 1. no efficient particle acceleration before $t_{pk}$
# 2. fixed fraction $\epsilon_{rel}$ of the shock power $L_{sh}$ is placed into relatistic particles
#
# Total energy of relativistic particles $E_{rel}$ is proportional to the fraction of radiated optical fluence $f_{sh}$:
#
# $
# E_{rel} \approx f_{sh} \epsilon_{rel} E_{opt}
# $
#
# $\epsilon_{rel} E_{opt}$ is an upper limit to the energy of accelerated particles, as $f_{sh}<1$.
#
# The total optical energy of shock-powered transients places an upper bound on the total gamma-ray and neutrino emission, as the relativistic particles cool quickly (**calorimetric method**).
#
# ## Coooling time scales
#
# (skip discussions on cooling time scales, **revisit later**)
#
# Shocks are generally radiative ($t_{cool} \ll t_{dyn}$) for $v_{sh} \lesssim 10,000-30,000$ km s$^{-1}$.
#
#
#
# For Novae: $\epsilon_{rel} \sim 0.003 - 0.01$
# **Shock dynamical time scale:**
#
# $t_{dyn} \sim R_{sh} / v_{sh}$
#
# Assume in the following $t_{dyn} \sim t_{pk}$.
#
# Approximation of opacity in fully ionized gas (electron scattering):
#
# $k_{es} \approx \sigma_T / m_p \approx 0.38 $ cm$^2$ g$^{-1}$.
#
# This is a good approximation for hydrogen-rich ejected; for hydro-poor ejecta, the opacity should be lower.
def t_dyn(n, v_sh, kappa = 0.38 * u.cm**2/u.g):
"dynamical time scale (peak time)"
return c/v_sh**2/n/(const.m_p*kappa)
# **Cooling times scale for radiatively cooled thermal gas**:
#
# $
# t_{cool} = \frac{9}{128} \frac{m_p v_{sh}^2}{\Lambda n_{sh}}
# $
#
# at $T = T_{sh}$ and with cooling function $\Lambda$.
#
# At high temperatures $T \gtrsim 10^{7.3}$ K with dominating free-free cooling:
#
# $\Lambda \approx \Lambda_{ff} \approx 2.3 \times 10^{-27} (T_{sh}/$K$)^{1/2}$ erg cm$^3$ s$^{-1}$
#
# Line emission cooling contributing at $10^5 < T < 10^{7.3}$:
#
# $\Lambda_{line} \approx 1.1 \times 10^{-22} (T_{sh}/$K$)^{-0.7}$ erg cm$^3$ s$^{-1}$
#
# (see discussion in Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press, 2011)
def t_cool(n, v_sh, T_sh):
"thermal cooling time scale"
lambda_ff=2.3e-27*math.sqrt(T_sh.to(u.k))
return 9.128*const.m_p*v_sh**2/lambda_ff/n_sh
# **Cooling time scale PP interactions**:
#
# $
# t_{pp} \approx \left( n \sigma_{pp} c \right)^{-1}
# $
#
# ($\sigma_{pp} = 5 \times 10^{-26}$ cm$^2$ at 1 PeV)
def t_pp(n, sigma_pp=5.e-25*u.cm**2):
"cooling time scale pp interactions"
return 1./n/sigma_pp/const.c
# **Cooling due to photo pion production:**
#
# Threshold for photo pion production:
#
# $
# E_{p,th} \approx (\epsilon_{p\gamma ,th}/\epsilon_{opt}) m_p c^2 = 1.4\times 10^{16} (\epsilon_{opt}/$10 eV$)^{-1}$ eV
#
# with $\epsilon_{p\gamma ,th}=(m_{\pi}+m_{\pi}^2/m_p)c^2 \approx 150$ MeV
# **Cooling time scale summary**
#
# All ratios at $t_{pk}$.
#
# Ratios $t_{cool}/t_{dyn}$ and $t_{pp}/t_{dyn}$ show that both thermal and non-thermal particles cool efficienctly around $t_{pk}$ and that a calometric approach can be considered for shock-powered transients.
# +
def v_85(v_sh):
"convenient shock speed"
return v_sh/(30000. * u.km/u.s)
def kappa_03(kappa):
"convenient kappa"
return kappa / (0.3 *u.cm**2/u.g)
def t_cool_dyn(v_sh, kappa=0.38*u.cm**2/u.g):
"ratio cooling to dynamical time scale"
return 1.e-3 * kappa_03(kappa) * v_85(v_sh)**3
def t_pp_dyn(v_sh, kappa=0.38*u.cm**2/u.g):
"ratio pp cooling to dynamical time scale"
return 1.e-3 * kappa_03(kappa) * v_85(v_sh)**2
def t_pp_pgamma(v_sh,e_opt=1*u.eV,f_omega=1.):
"ratio of pp to p_gamma cooling"
return 4.*v_85(v_sh)**2/e_opt*f_omega
v_sh = (np.logspace( 3, 5.5, 100))*u.km/u.s
plt.gcf().clear()
plt.figure(figsize=(6,6))
r_cool_dyn = t_cool_dyn(v_sh)
plt.plot(
v_sh,
r_cool_dyn,
label="$t_{cool}/t_{dyn}$")
# H poor environement with lower opacity
r_cool_dyn_Hpoor = t_cool_dyn(v_sh, 0.5*0.38*u.cm**2/u.g)
#plt.plot(
# v_sh,
# r_cool_dyn_Hpoor,
# label="$t_{cool}/t_{dyn}$ (H poor)")
r_pp_dyn = t_pp_dyn(v_sh)
plt.plot(
v_sh,
r_pp_dyn,
label="$t_{pp}/t_{dyn}$")
# depends strongly on f_omega!
r_pp_pgamma = t_pp_pgamma(v_sh, 1., 1.)
plt.plot(
v_sh,
r_pp_pgamma,
label="$t_{pp}/t_{p\gamma}$")
plt.xlabel("shock velocity ($km/s$)",fontsize=14)
plt.ylabel("ratio of cooling time scales",fontsize=14)
plt.yscale('log')
plt.xscale('log')
plt.legend()
plt.title("Ratio of cooling time scales" )
plt.show()
# -
# ## Diffusive shock acceleration
#
# **Maximum Ion Energy:**
#
# Magnetic field strenght near shock (derived using equipartion arguments?):
#
# $
# B_{sh} = \left( 6 \pi \epsilon_B m_p n_{sh} v_{sh}^2 \right)^{1/2}
# $
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment