Source code for gammapy.astro.source.pulsar

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Pulsar source models."""
import numpy as np
from astropy.units import Quantity
from ...extern.validator import validate_physical_type

__all__ = ["Pulsar", "SimplePulsar"]

DEFAULT_I = Quantity(1e45, "g cm2")
"""Pulsar default moment of inertia"""

DEFAULT_R = Quantity(1e6, "cm")
"""Pulsar default radius of the neutron star"""

B_CONST = Quantity(3.2e19, "gauss s^(-1/2)")
"""Pulsar default magnetic field constant"""


[docs]class SimplePulsar: """Magnetic dipole spin-down model for a pulsar. Reference: http://www.cv.nrao.edu/course/astr534/Pulsars.html Parameters ---------- P : `~astropy.units.Quantity` Rotation period (sec) P_dot : `~astropy.units.Quantity` Rotation period derivative (sec sec^-1) I : `~astropy.units.Quantity` Moment of inertia (g cm^2) R : `~astropy.units.Quantity` Radius of the pulsar (cm) """ def __init__(self, P, P_dot, I=DEFAULT_I, R=DEFAULT_R): validate_physical_type("P", P, "time") validate_physical_type("P_dot", P_dot, "dimensionless") self.P = P self.P_dot = P_dot self.I = I self.R = R @property def luminosity_spindown(self): """Spin-down luminosity (`~astropy.units.Quantity`). Notes ----- The spin-down luminosity is given by: .. math:: \\dot{L} = 4\\pi^2 I \\frac{\\dot{P}}{P^{3}} """ return 4 * np.pi ** 2 * self.I * self.P_dot / self.P ** 3 @property def tau(self): """Characteristic age (`~astropy.units.Quantity`). Notes ----- The characteristic age is given by: .. math:: \\tau = \\frac{P}{2\\dot{P}} """ return (self.P / (2 * self.P_dot)).to("yr") @property def magnetic_field(self): """Magnetic field strength at the polar cap (`~astropy.units.Quantity`). Notes ----- The magnetic field is given by: .. math:: B = 3.2\\cdot 10^{19} (P\\dot{P})^{1/2} [\\textnormal(Gauss)] """ return B_CONST * np.sqrt(self.P * self.P_dot)
[docs]class Pulsar(SimplePulsar): """Magnetic dipole spin-down pulsar model. Reference: http://www.cv.nrao.edu/course/astr534/Pulsars.html Parameters ---------- P_0 : float Period at birth logB : float Logarithm of the magnetic field, which is constant n : float Spin-down braking index I : float Moment of inertia R : float Radius """ def __init__( self, P_0=Quantity(0.1, "s"), logB=10, n=3, I=DEFAULT_I, R=DEFAULT_R, age=None, L_0=None, morphology="Delta2D", ): self.I = I self.R = R self.P_0 = P_0 self.logB = logB self.P_dot_0 = (Quantity(10 ** logB, "gauss") / B_CONST) ** 2 / P_0 self.tau_0 = P_0 / (2 * self.P_dot_0) self.n = float(n) self.beta = (n + 1.0) / (n - 1.0) self.morphology = morphology if age is not None: validate_physical_type("age", age, "time") self.age = age if L_0 is None: self.L_0 = 4 * np.pi ** 2 * self.I * self.P_dot_0 / self.P_0 ** 3
[docs] def luminosity_tev(self, t=None, fraction=0.1): """Gamma-ray luminosity assumed to be a certain fraction of the spin-down luminosity. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return self.luminosity_spindown(t) * fraction
[docs] def luminosity_spindown(self, t=None): """Spin down luminosity at age t. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. Notes ----- The spin-down luminosity is given by: .. math:: \\dot{L}(t) = \\dot{L}_0 \\left(1 + \\frac{t}{\\tau_0}\\right)^{\\frac{n + 1}{n - 1}} """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return self.L_0 * (1 + (t / self.tau_0)) ** self.beta
[docs] def period(self, t=None): """Period at age t. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. Notes ----- The period is given by: .. math:: P(t) = P_0\\left(1 + \\frac{t}{\\tau_0}\\right)^{\\frac{1}{n - 1}} """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return self.P_0 * (1 + (t / self.tau_0)) ** self.beta
[docs] def energy_integrated(self, t=None): """Total released energy at age t. Time-integrated spin-down luminosity since birth. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. Notes ----- The time integrated energy is given by: .. math:: E(t) = \\dot{L}_0 \\tau_0 \\frac{t}{t + \\tau_0} """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return self.L_0 * self.tau_0 * (t / (t + self.tau_0))
[docs] def period_dot(self, t=None): """Period derivative at age t. P_dot for a given period and magnetic field B, assuming a dipole spin-down. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. Notes ----- The period derivative is given by: .. math:: \\dot{P}(t) = \\frac{B^2}{3.2 \\cdot 10^{19} P(t)} """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return Quantity(10 ** self.logB, "gauss") ** 2 / (self.period(t) * B_CONST ** 2)
[docs] def tau(self, t=None): """Characteristic age at real age t. Parameters ---------- t : `~astropy.units.Quantity` Time after birth of the pulsar. Notes ----- The characteristic age is given by: .. math:: \\tau = \\frac{P}{2\\dot{P}} """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return self.period(t) / 2 * self.period_dot(t)
[docs] def magnetic_field(self, t=None): """Magnetic field strength at the polar cap. Assumed to be constant. Notes ----- The magnetic field is given by: .. math:: B = 3.2\\cdot 10^{19} (P\\dot{P})^{1/2} [\\textnormal(Gauss)] """ if t is not None: validate_physical_type("t", t, "time") elif hasattr(self, "age"): t = self.age else: raise ValueError("Need time variable or age attribute.") return B_CONST * np.sqrt(self.period(t) * self.period_dot(t))