Spikes in IC spectra at highest energies

[alziegler]

Hi @zblz ,

thank you for providing this nice modeling tool!

I used the code from your documentation (naima version 0.8), to reproduce the IC plot shown here: http://naima.readthedocs.io/en/latest//radiative-1.png

I modified the range shown on the y-axis and encountered some irregularities in the NIR IC spectra at the highest energies, which you can see in the following plots:

After that, I changed the energy range up to 1e18 and switched from cutoff power law to pure power law model for the electron spectrum, yielding the following spectra:

I wonder where the spikes in the spectra at the highest energies are coming from? Am I doing something wrong in the calling of the functions/computation (but I did not really change something at all)? Or is it related to the Klein Nishina supression, which one after another of the seed fields encounter at higher energies (related to the spectral breaks one can observe), which could eventually point to effects due to simplifications/assumptions made in the computation of the spectra?

In the following the code, used to produce the plots:

import naima
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt

# Define models
#spectral_model = naima.models.ExponentialCutoffPowerLaw(1e36*u.Unit('1/eV'),
#        1*u.TeV, 2.1, 13*u.TeV)

spectral_model = naima.models.PowerLaw(1e36*u.Unit('1/eV'), 1*u.TeV, 2.0)

IC = naima.models.InverseCompton(spectral_model, seed_photon_fields=['CMB', 'FIR', 'NIR'], Eemin=1*u.GeV, Eemax=1.e18*u.eV)
#IC.particle_distribution.index = 1.8
SYN = naima.models.Synchrotron(spectral_model, B=100*u.uG)

# Compute SEDs
spectrum_energy = np.logspace(-1,18,1000)*u.eV
sed_IC = IC.sed(spectrum_energy, distance=1.5*u.kpc)
sed_SYN = SYN.sed(spectrum_energy, distance=1.5*u.kpc)

# Plot
plt.figure(figsize=(8,5))
plt.rc('font', family='sans')
plt.rc('mathtext', fontset='custom')
for seed, color in zip(['CMB', 'FIR', 'NIR'], ['lime','black','red']):
    sed = IC.sed(spectrum_energy, seed=seed, distance=1.5*u.kpc)
    plt.loglog(spectrum_energy,sed,lw=1,
               ls='-', c=color, label='IC ({0})'.format(seed))
plt.loglog(spectrum_energy,sed_IC,lw=2,
        label='IC (total)',c=naima.plot.color_cycle[0])
plt.loglog(spectrum_energy,sed_SYN,lw=2,label='Sync',c=naima.plot.color_cycle[1])
plt.xlabel('Photon energy [{0}]'.format(
        spectrum_energy.unit.to_string('latex_inline')))
plt.ylabel('$E^2 dN/dE$ [{0}]'.format(
        sed_SYN.unit.to_string('latex_inline')))
plt.ylim(1e-18, 1e-6)
plt.tight_layout()
plt.legend(loc='lower left')

plt.show()

[zblz]

Hi @alziegler - This is an expected side effect of the way IC is computed. The particle distribution is sampled in log-space and the IC contribution from each energy computed and integrated numerically. When IC is in deep KN, the emission profile for a given electron energy is very sharp: effectively all of the electron energy is transfered to the photon, so the spectrum is quasi-monochromatic. The combination of the way the spectrum is computed and this effect in deep KN results in the spikes you see.

These can be avoided by sampling more densely the electron spectrum, but this also results in a longer computation time, so the defaults are set quite low. You can increase the density using the nEed parameter of the InverseCompton class from its default of 300 points per energy decade (see the Other Parameters section of the API docs).

[alziegler]

@zblz - thank you for the quick answer. I switched to larger nEed values, which, as you explain, improves the calculation artifacts. It might also be helpful for other users to hint at this potential problem somewhere in the docs? So I think the issue can be closed right away.