Cutoff in the SSC spectrum at high frequencies

[cosimoNigro]

Dear @andreatramacere,

I was making some comparison plots between agnpy, jetset (version 1.1.2) and some reference data I had. The SSC of jetset always shows a cutoff at high frequencies, compared to agnpy and to the reference SED.

Can you clarify to what it is due and how can I obtain a smoother SSC SED? Here the snippet of code I use to compute that SED - parameters are set accordingly to the reference SED:

# generate a synchrotron + SSC SED to be confronted with the one produced by
# agnpy and Figure 7.4 of Dermer 2009
from jetset.jet_model import Jet
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt

# jet with power-law electron distribution
pwl_jet = Jet(
    name="", electron_distribution="pl", electron_distribution_log_values=False
)
# set parameters according to Figure 7.4 of Dermer 2009
pwl_jet.set_par("N", 1298.13238394)
pwl_jet.set_par("p", 2.8)
pwl_jet.set_par("gmin", 1e2)
pwl_jet.set_par("gmax", 1e5)
pwl_jet.set_par("B", 1)
pwl_jet.set_par("R", 1e16)
pwl_jet.set_par("beam_obj", 10)
pwl_jet.set_par("z_cosm", 0.07)
# remove SSA
pwl_jet.spectral_components.Sync.state = "on"
# SSC emission
pwl_jet.set_nu_grid(1e14, 1e26, 100)
pwl_jet.show_model()
pwl_jet.eval()

ssc_nu = pwl_jet.spectral_components.SSC.SED.nu
ssc_sed = pwl_jet.spectral_components.SSC.SED.nuFnu
plt.loglog(ssc_nu, ssc_sed)
plt.ylim([1e-20, 1e-9])
plt.show()

is there some cutoff you are imposing in the numerical calculation?

As always, thanks a lot!

P.S. maybe "cutoff" is not the proper term, I just mean that the SED flux goes abruptly to zero after a given frequency.

[andreatramacere]

Hi, you should increase the number of points in the IC component, in this version the default version is 50. I am attaching a script done with the pre-release of 1.2.0 that mimics 1.1.2, and that should work with 1.1.2. Jetset reproduces all the spectral components on the same grid (for a given model), and points below a given source emissivity lower bound are set to the lower bound value (the default is 1E-120). Having the same grids allows to combine the spectral components easily. This explain the sharp 'cut-off' in the plot. In v 1.2.0 these settings are more relaxed for the user. By the way, you should test with the new version that I will release soon "v.1.2.0". You can send me the tests and I will send you back the results fo this and other scripts for v.1.2.0. Or if you want you can use the pre-release. Also checking the actual synchrotron kernel might be relevant to these tests.

# generate a synchrotron + SSC SED to be confronted with the one produced by
# agnpy and Figure 7.4 of Dermer 2009
from jetset.jet_model import Jet
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt

# jet with power-law electron distribution
pwl_jet = Jet(
    name="", electron_distribution="pl", electron_distribution_log_values=False
)
# set parameters according to Figure 7.4 of Dermer 2009
pwl_jet.set_par("N", 1298.13238394)
pwl_jet.set_par("p", 2.8)
pwl_jet.set_par("gmin", 1e2)
pwl_jet.set_par("gmax", 1e5)
pwl_jet.set_par("B", 1)
pwl_jet.set_par("R", 1e16)
pwl_jet.set_par("beam_obj", 10)
pwl_jet.set_par("z_cosm", 0.07)
# remove SSA
pwl_jet.spectral_components.Sync.state = "on"
# SSC emission
pwl_jet.set_nu_grid(1e14, 1e26, 100)
#to mimic old default behaviuor in 1.1.2
pwl_jet.set_IC_nu_size(50)
pwl_jet.show_model()
pwl_jet.eval()

ssc_nu = pwl_jet.spectral_components.SSC.SED.nu
ssc_sed = pwl_jet.spectral_components.SSC.SED.nuFnu
plt.loglog(ssc_nu, ssc_sed,label='50 IC points')
plt.ylim([1e-13, 1e-9])
#new default behaviuor in 1.2.0
pwl_jet.set_IC_nu_size(100)
pwl_jet.eval()
ssc_nu = pwl_jet.spectral_components.SSC.SED.nu
ssc_sed = pwl_jet.spectral_components.SSC.SED.nuFnu
plt.loglog(ssc_nu, ssc_sed,label='100 IC points')
#to remove points below flux_plot_lim
msk= pwl_jet.spectral_components.SSC.SED.nuFnu.value>pwl_jet.flux_plot_lim
ssc_nu = pwl_jet.spectral_components.SSC.SED.nu
ssc_sed = pwl_jet.spectral_components.SSC.SED.nuFnu
plt.loglog(ssc_nu[msk], ssc_sed[msk],label='100 IC points, flux>flux_plot_lim')
plt.legend()
plt.show()

[cosimoNigro]

Thanks for the clarification @andreatramacere,

I regenerated the plot with a larger IC_nu_grid.

I have a script that I run manually to generate the jetset reference SEDs. https://github.com/cosimoNigro/agnpy/blob/master/agnpy/data/reference_seds/jetset/jetset_ssc_seds.py The tests of agnpy take care of comparing the two automatically.

So I will re-generate the jetset reference SEDs once you make the 1.2.0 release. I am not in a hurry.

You can close this, thanks.