Naima Radiative Models

Welcome to the naima radiative models tutorial!

Useful references:

2nd part of this tutorial: MCMC model fitting of a galactic non-thermal source

naima is named after a ballad composed by John Coltrane, so a few performances by him and others might be an appropriate soundtrack for the tutorial.



In [ ]:

    
#prepare imports
import numpy as np
import astropy.units as u
import naima
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['lines.linewidth'] = 2

Note: Astropy units are used throughout this tutorial, see a short primer in the naima docs and a longer one in the astropy docs.

Exploring the naima radiative models

naima provides several radiative models and functional models than can be used together to compute the radiative output from populations of relativistic particle populations.

General concepts

In general, the radiative output for a given channel can be computed as:

$$L(\epsilon)=\iiiint E\,N(E, \vec{r})\,c\,n(\epsilon_0)\,\frac{\mathrm{d} \sigma(E,\epsilon_0,\epsilon,\Omega)}{\mathrm{d} \Omega} \mathrm{d}\Omega\, \mathrm{d}\epsilon_0\, \mathrm{d} E\, \mathrm{d} V, $$

where $E$ is the particle energy, $\epsilon$ the emitted photon energy, $\sigma$ the cross-section of the process, which might depend on the angle of interaction, $n(\epsilon_0)$ is the density of target particles or photons at energy $\epsilon_0$, and $N(E)$ is the energy distribution of the relativistic particle population. In a one-zone model, the properties are assumed homogeneous over the whole volume, so one fewer integral is needed:

$$L(\epsilon)=\iiint E\,N(E)\,c\,n(\epsilon_0)\,\frac{\mathrm{d} \sigma(E,\epsilon_0,\epsilon,\Omega)}{\mathrm{d} \Omega} \mathrm{d}\Omega\, \mathrm{d}\epsilon_0\, \mathrm{d} E, $$

One of the most prevalent mechanisms of particle acceleration is diffusive shock acceleration, which, for a strong shock, results in a powerlaw particle energy distribution with index $p=2$. The acceleration timescale increases with particle energy, and at some energy it will be longer than either the radiative or non-radiative losses of the particle or longer than the age of the source. At this energy the particle distribution will show a quasi-exponential cutoff.

Energy losses might modify the present-age particle distribution, typically making the distribution softer. This might only happen at higher energies, giving rise to a broken power-law energy distribution.

For a given relativistic particle population, there are different channels that can result in photon emission depending on the target of the interaction:

Synchrotron: Charged particles will radiate as they girate in magnetic fields.
Inverse Compton: Upscattering of seed photon fields by electrons.
Bremsstrahlung: Charged particles radiate as they are accelerated by nearby particles.
Pion Decay: Proton-proton interactions results in pions that decay into gamma-rays.

naima models

naima has the following radiative models available in the naima.models module:

Synchrotron
InverseCompton
PionDecay
Bremsstrahlung

All of them take as first argument a particle distribution function, which can be one of the functional models in naima (currently PowerLaw, BrokenPowerLaw, ExponentialCutoffPowerLaw, ExponentialCutoffBrokenPowerLaw, LogParabola). You can find out the arguments they take by looking at their docstrings:



In [ ]:

    
naima.models.Synchrotron?

Leptonic models

The first step is to define the particle distribution we will use, and we'll start with a power-law with an exponential cutoff. The amplitude must be in units of particles per unit energy:



In [ ]:

    
ECPL = naima.models.ExponentialCutoffBrokenPowerLaw(amplitude=1e36/u.eV, e_0=1*u.TeV, alpha_1=2.0,
                                                    alpha_2=3.0, e_break=0.2*u.TeV, e_cutoff=50*u.TeV,
                                                    beta=2)
#check the shape of the particle distribution
electron_energy = np.logspace(-3, 3, 1000) * u.TeV
plt.loglog(electron_energy, (electron_energy**2 * ECPL(electron_energy)).to('erg'))
plt.axvline(0.2, ls=':', lw=1, c='k'); plt.axvline(50, ls=':', lw=1, c='k')
plt.gca().set_ylim(bottom=1e42)
plt.xlabel('Particle energy [TeV]')
plt.ylabel('$E^2 N(E)$ [erg]')

Now we can define the radiative channels we want to compute with this particle distribution. We will consider it to be an electron population, and compute its synchrotron SED in a magnetic field strength of $B=10\mu G$ at a distance of 1 kpc:



In [ ]:

    
SY = naima.models.Synchrotron(ECPL, B=10*u.uG)
photon_energy = np.logspace(-5, 5, 100)*u.eV
SY_sed0 = SY.sed(photon_energy, distance=1*u.kpc)
plt.loglog(photon_energy, SY_sed0)
plt.axvspan(1e3,1e4, fc='k', alpha=0.3)
plt.xlabel('Photon energy [eV]')
plt.ylabel('$E^2 dN/dE\ [\mathrm{erg\ cm^{-2}\ s^{-1}}]$')

Now we can use the same particle distribution in an Inverse Compton class to compute its IC emission. We will use the interstellar radiation fields computed from the GALPROP model at a galactocentric distance of ~6 kpc:

Name	Temperature	Energy density
CMB	2.72 K	0.216 eV/cm$^3$
Far infrared	27 K	0.415 eV/cm$^3$
Near infrared	2800 K	0.802 eV/cm$^3$

The total IC emission should be slighlty lower than the synchrotron emission, given that (in the Thomson regime for IC) the ratio between the Synchrotron and IC luminosities corresponds to the ration of target energy densities:

$$\frac{L_{sy}}{L_{IC}} = \frac{u_B}{u_\gamma} = \frac{B^2}{8\pi u_\gamma} \simeq 1.73$$



In [ ]:

    
IC = naima.models.InverseCompton(ECPL,
                   seed_photon_fields=['CMB',
                                       ['FIR', 27*u.K, 0.415*u.eV/u.cm**3],
                                       ['NIR', 2800*u.K, 0.802*u.eV/u.cm**3]])

# The particle_distribution attribute in IC and SY now point to the same object:
print('Are the SY and IC particle distributions the same object? {0}'.format(
        IC.particle_distribution is SY.particle_distribution))

# and we compute the SED from optical to TeV
photon_energy = np.logspace(-5, 15, 100) * u.eV
f, ax = plt.subplots(1)

# plot the total SED from SY and IC
ax.loglog(photon_energy, SY.sed(photon_energy), label='Synchrotron')
ax.loglog(photon_energy, IC.sed(photon_energy), label='IC (total)')

# plot the SEDs from each of the seed photon fields with seed argument
ax.loglog(photon_energy, IC.sed(photon_energy, seed='CMB'), 
              label='IC (CMB)', ls='--', lw=1)
ax.loglog(photon_energy, IC.sed(photon_energy, seed='FIR'), 
              label='IC (FIR)', ls='--', lw=1)
ax.loglog(photon_energy, IC.sed(photon_energy, seed='NIR'), 
              label='IC (NIR)', ls='--', lw=1)

ax.set_xlabel('Photon energy [eV]')
ax.set_ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')
ax.legend(loc='lower left')
ax.set_ylim(bottom=1e-14)

Now we will explore how to modify the parameters of an already defined radiative model. We can modify the values of the particle distribution parameters and this will propagate to the radiative classes:



In [ ]:

    
ECPL.e_cutoff = 50*u.TeV
plt.loglog(photon_energy, SY.sed(photon_energy) + IC.sed(photon_energy), label='$E_\mathrm{c}$ = 50 TeV')
ECPL.e_cutoff = 100*u.TeV
plt.loglog(photon_energy, SY.sed(photon_energy) + IC.sed(photon_energy), label='$E_\mathrm{c}$ = 100 TeV')

plt.legend(loc='best')
plt.xlabel('Photon energy [eV]')
plt.ylabel('$E^2 dN/dE$ [erg/cm2/s]')
plt.ylim(bottom=1e-12)

You can also modify the normalization by specifying the total energy in electrons with the set_We function:



In [ ]:

    
f, ax = plt.subplots(1)
print('before setting We -- ECPL.amplitude = {0}'.format(ECPL.amplitude))
IC.set_We(5e48 * u.erg, Eemin=0.1*u.TeV) # 5x10^48 erg in electrons above 100 GeV
print('after  setting We -- ECPL.amplitude = {0}'.format(ECPL.amplitude))
ax.loglog(photon_energy, SY.sed(photon_energy) + IC.sed(photon_energy), label=r'$5\times10^{48}$ erg', c='k')

IC.set_We(5e47 * u.erg, Eemin=0.1*u.TeV) # 5x10^47 erg in electrons above 100 GeV
ax.loglog(photon_energy, SY.sed(photon_energy) + IC.sed(photon_energy), label=r'$5\times10^{47}$ erg', c='r')

#plot label stuff
ax.set_xlabel('Photon energy [eV]')
ax.set_ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')
ax.legend(loc='lower left')
ax.set_ylim(bottom=1e-13)



In [ ]:

    
f, ax = plt.subplots(1)
for B in [0.1, 1, 10, 100, 1000]*u.uG:
    SY.B = B
    ax.loglog(photon_energy, SY.sed(photon_energy) + IC.sed(photon_energy),
              label=r'$B$ = {0}'.format(B))
    
#plot label stuff
ax.set_xlabel('Photon energy [eV]')
ax.set_ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')
ax.legend(loc='best')
ax.set_ylim(bottom=1e-15)

Klein-Nishina demo

The regime of Inverse Compton radiation depends on the product of the energy of the seed photon $\epsilon_0$ and the energy of the electron $E_e$ in units of $m_e c^2$: $\kappa_0=\epsilon_0 E_e/(m_e^2 c^4)$. For $\kappa_0<<1$ IC proceeds in the Thomson regime (with cross section increasing with energy), and at $\kappa_0>>1$ it proceeds in the Klein-Nishina regime, in which the emission from a gievn electron energy is close to monochromatic and the cross section diminishes rapidly with energy, leading to a break in the photon spectrum in the transition between the two regimes.



In [ ]:

    
PL = naima.models.PowerLaw(1e36/u.eV, e_0=1*u.TeV, alpha=3.0)
IC_KN = naima.models.InverseCompton(PL, seed_photon_fields=[
                                        ['CMB', 2.7*u.K,  1*u.eV/u.cm**3],
                                        ['FIR', 70*u.K,   1*u.eV/u.cm**3],
                                        ['NIR', 1800*u.K, 1*u.eV/u.cm**3]])

# Set maximum and minimum electron energies to a very wide range
IC_KN.Eemax = 10*u.PeV
IC_KN.Eemin = 1*u.MeV

# compute the three IC spectra
photon_energy = np.logspace(-6, 3, 100) * u.TeV
for seed in ['CMB', 'FIR', 'NIR']:
    plt.loglog(photon_energy, IC_KN.sed(photon_energy, seed=seed), 
              label='$T_{{0}}$ = {0:.1f}'.format(IC.seed_photon_fields[seed]['T']), lw=1)
    
# plot labels
plt.legend(loc='best')
plt.ylim(bottom=1e-12)
plt.xlabel('Photon energy [TeV]')
plt.ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')

Pion Decay

Now we can explore the pion decay emission from a relativistic proton distribution. We define a new particle distribution:



In [ ]:

    
proton_dist = naima.models.ExponentialCutoffPowerLaw(1e36/u.eV, e_0=1*u.TeV, alpha=2.1, e_cutoff=100*u.TeV)
PP = naima.models.PionDecay(proton_dist, nh=0.1*u.cm**-3)

photon_energy = np.logspace(-2, 6, 100) * u.GeV
plt.loglog(photon_energy, PP.sed(photon_energy))

#plot label stuff
plt.xlabel('Photon energy [GeV]')
plt.ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')
plt.ylim(bottom=1e-14)

You will need to install ipywidgets for this to work



In [ ]:

    
ECPL = naima.models.ExponentialCutoffPowerLaw(amplitude=1e36/u.eV, e_0=1*u.TeV, alpha=2.5,
                                                    e_cutoff=50*u.TeV, beta=2)

IC = naima.models.InverseCompton(ECPL, seed_photon_fields=['CMB', ['FIR', 90*u.K, 0.415*u.eV/u.cm**3]])
#IC.nEed /=10
SY = naima.models.Synchrotron(ECPL, B=10*u.uG)
photon_energy = np.logspace(-5, 15, 100)*u.eV

from ipywidgets import interact
from IPython.display import display

@interact(logA=(33.,39.,0.2), alpha=(1.,4.,0.1), log_e_cutoff=(-1,3,0.2), u_FIR=(0.1,5.,0.1),  T_FIR=(10,300,10), logB=(-1,2,0.1))
def model(logA=36, alpha=2.5, log_e_cutoff=1.8, u_FIR=0.4, T_FIR=90, logB=1):
    ECPL.amplitude = 10**logA /u.eV
    ECPL.alpha = alpha
    ECPL.e_cutoff = 10**log_e_cutoff * u.TeV
    IC.seed_photon_fields['FIR']['u'] = u_FIR * u.eV/u.cm**3
    IC.seed_photon_fields['FIR']['T'] = T_FIR * u.K
    SY.B = 10**logB*u.uG
    
    f,ax = plt.subplots(1)
    ax.loglog(photon_energy, SY.sed(photon_energy))
    ax.loglog(photon_energy, IC.sed(photon_energy))
    ax.loglog(photon_energy, IC.sed(photon_energy, seed='CMB'), lw=1, ls='--')
    ax.loglog(photon_energy, IC.sed(photon_energy, seed='FIR'), lw=1, ls='--')
    plt.axvspan(1e2,1e4, fc='0.5', alpha=0.25, lw=0)
    plt.axvspan(1e8,1e10, fc='r', alpha=0.25, lw=0)
    plt.axvspan(1e11,1e13, fc='b', alpha=0.25, lw=0)
    ax.set_ylim(bottom=1e-13)
    ax.set_xlabel('Photon energy [eV]')
    ax.set_ylabel('$E^2 dN/dE$ [erg s$^{-1}$ cm$^{-2}$]')
    plt.show()
    print('W_e(E_e>1 TeV) = {0:.3g}'.format(IC.compute_We(Eemin=1*u.TeV)))
    print('E_c = {0:.3g}'.format(ECPL.e_cutoff))
    print('B = {0:.3g}'.format(SY.B))
    print('u_FIR = {0:.3g}'.format(IC.seed_photon_fields['FIR']['u']))
    print('T_FIR = {0:.3g}'.format(IC.seed_photon_fields['FIR']['T']))