# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Pulsar wind nebula (PWN) source models."""
import numpy as np
from astropy.units import Quantity
from astropy.utils import lazyproperty
import astropy.constants as const
from scipy.optimize import fsolve
from ..source import Pulsar, SNRTrueloveMcKee
__all__ = ["PWN"]
[docs]class PWN:
"""Simple pulsar wind nebula (PWN) evolution model.
Parameters
----------
pulsar : `~gammapy.astro.source.Pulsar`
Pulsar model instance.
snr : `~gammapy.astro.source.SNRTrueloveMcKee`
SNR model instance
eta_e : float
Fraction of energy going into electrons.
eta_B : float
Fraction of energy going into magnetic fields.
age : `~astropy.units.Quantity`
Age of the PWN.
morphology : str
Morphology model of the PWN
"""
def __init__(
self,
pulsar=Pulsar(),
snr=SNRTrueloveMcKee(),
eta_e=0.999,
eta_B=0.001,
morphology="Gaussian2D",
age=None,
):
self.pulsar = pulsar
if not isinstance(snr, SNRTrueloveMcKee):
raise ValueError("SNR must be instance of SNRTrueloveMcKee")
self.snr = snr
self.eta_e = eta_e
self.eta_B = eta_B
self.morphology = morphology
if age is not None:
self.age = Quantity(age, "yr")
def _radius_free_expansion(self, t):
"""Radius at age t during free expansion phase.
Reference: http://adsabs.harvard.edu/abs/2006ARA%26A..44...17G (Formula 8).
"""
term1 = (self.snr.e_sn ** 3 * self.pulsar.L_0 ** 2) / (self.snr.m_ejecta ** 5)
return (1.44 * term1 ** (1.0 / 10) * t ** (6.0 / 5)).cgs
@lazyproperty
def _collision_time(self):
"""Time of collision between the PWN and the reverse shock of the SNR.
Returns
-------
t_coll : `~astropy.units.Quantity`
Time of collision.
"""
def time_coll(t):
t = Quantity(t, "yr")
r_pwn = self._radius_free_expansion(t).to_value("cm")
r_shock = self.snr.radius_reverse_shock(t).to_value("cm")
return r_pwn - r_shock
# 4e3 years is a typical value that works for fsolve
return Quantity(fsolve(time_coll, 4e3), "yr")
[docs] def radius(self, t):
"""Radius of the PWN at age t.
During the free expansion phase the radius of the PWN evolves like:
.. math::
R_{PWN}(t) = 1.44 \\left(\\frac{E_{SN}^3\\dot{E}_0^2}
{M_{ej}^5}\\right)^{1/10}t^{6/5}
\\text{pc}
After the collision with the reverse shock of the SNR, the radius is
assumed to be constant (See `~gammapy.astro.source.SNRTrueloveMcKee.radius_reverse_shock`).
Reference: http://adsabs.harvard.edu/abs/2006ARA%26A..44...17G (Formula 8).
Parameters
----------
t : `~astropy.units.Quantity`
Time after birth of the SNR
"""
t = Quantity(t, "yr")
r_collision = self._radius_free_expansion(self._collision_time)
r = np.where(
t < self._collision_time,
self._radius_free_expansion(t).value,
r_collision.value,
)
return Quantity(r, "cm")
[docs] def magnetic_field(self, t):
"""Estimate of the magnetic field inside the PWN.
By assuming that a certain fraction of the spin down energy is
converted to magnetic field energy an estimation of the magnetic
field can be derived.
Parameters
----------
t : `~astropy.units.Quantity`
Time after birth of the SNR
"""
t = Quantity(t, "yr")
energy = self.pulsar.energy_integrated(t)
volume = 4.0 / 3 * np.pi * self.radius(t) ** 3
return np.sqrt(2 * const.mu0 * self.eta_B * energy / volume)