%load_ext autoreload
%autoreload 2

from astropy import units
import numpy as np

Flux Definition

SkyLLH provides a sophisticated class collection to describe a differential particle flux function.

General Flux Function

The most general differential flux function is:

\(\frac{\mathrm{d}^4\Phi(\alpha,\delta,E,t | \vec{p}_\mathrm{s})}{\mathrm{d}A \mathrm{d}\Omega \mathrm{d}E \mathrm{d}t}\),

which is a function of celestial coordinates right-ascention, \(\alpha\), and declination, \(\delta\), energy \(E\), and time \(t\), given source parameters \(\vec{p}_\mathrm{s}\), e.g. source location and spectral index \(\gamma\) for a power-law energy profile.

The abstract base class for all flux function models is FluxModel, which is derived from MathFunction and Model.

The FluxModel has the abstract __call__() method defined, which will evaluate the flux function for values of \(\alpha\), \(\delta\), \(E\), and \(t\) given in specific units.

Units for flux models

The FluxModel class defines the units used for length, angle, energy, and time. The default units are configured through a skyllh.core.config.Config instance, in particular through the cfg['units']['defaults']['fluxes'] dictionary of the Config instance. Units must be derived classes from the astropy.units.UnitBase class.

The FluxModel class has the properties unit_str and unit_latex_str for a representation of the units as a str object in plain text and latex code.

Factorized Flux Function

Usually the flux function can be split into a spatial, energy, and time profile with an overall normalization constant in differential flux unit, i.e. \(\mathrm{area}^{-1} \mathrm{solid-angle}^{-1} \mathrm{energy}^{-1} \mathrm{time}^{-1}\). Hence, SkyLLH provides the class FactorizedFluxModel, which describes a differential flux function as the product of individual flux profiles:

\[\Phi(\alpha,\delta,E,t | \vec{p}_\mathrm{s}) = \Phi_0 \Psi(\alpha,\delta|\vec{p}_\mathrm{s}) \epsilon(E|\vec{p}_\mathrm{s}) T(t|\vec{p}_\mathrm{s})\]

The abstract base class for any flux profile is FluxProfile, which is derived from MathFunction.

The abstract base class for a spatial, energy, and time flux profile is SpatialFluxProfile, EnergyFluxProfile, and TimeFluxProfile, respectively, and are derived from FluxProfile.

Steady Point-Like Flux

A very common source hypothesis is a steadily emitting point-like source. Hence, SkyLLH provides the class SteadyPointlikeFFM. It takes a flux normalization \(\Phi_0\), and an energy profile as constructor arguments. As spatial profile it uses the PointSpatialFluxProfile class.

As an example we create a steady point-like factorized flux model with a power-law energy flux profile.

But first we need to create a configuration:

from skyllh.core.config import Config
cfg=Config()

Now we see what the default units are:

from skyllh.core.flux_model import FluxModel

print(FluxModel.get_default_units(cfg=cfg))

{'angle': Unit("rad"), 'energy': Unit("GeV"), 'length': Unit("cm"), 'time': Unit("s")}

Now we need to create the energy flux profile. As reference energy and spectral index we choose $E_0=10^3$~GeV and $\gamma=2$, respectively:

from skyllh.core.flux_model import PowerLawEnergyFluxProfile

energy_profile = PowerLawEnergyFluxProfile(
    E0=1e3,
    gamma=2,
    energy_unit=units.GeV,
    cfg=cfg)
print(energy_profile)

(E / (1000 GeV))^-2

The next step is to create the SteadyPointlikeFFM class instance. As normalization constant we choose \(\Phi_0 = 10^{-8} \text{GeV}^{-1}\text{cm}^{-2}\text{sr}^{-1}\text{s}^{-1}\):

from skyllh.core.flux_model import SteadyPointlikeFFM

fluxmodel = SteadyPointlikeFFM(
    Phi0=1e-8,
    energy_profile=energy_profile,
    angle_unit=units.radian,
    time_unit=units.s,
    length_unit=units.cm,
    cfg=cfg
)
print(fluxmodel)

1.000e-08 * (E / (1000 GeV))^-2 * 1 (GeV cm^2 s)^-1

Evaluating the flux model function

The flux model function can be evaluated by calling its __call__() operator:

flux = fluxmodel(E=3e3)
print(f'flux.shape = {flux.shape}')
print(f'flux = {flux}')

flux.shape = (1, 1, 1)
flux = [[[1.11111111e-09]]]

It returns a 3-dimensional numpy array, where the first, second, and third dimension represents the spatial, energy, and time axes, respectively. Hence, we can evaluate the flux model for different spatial coordinates, energies, and times by a single call:

flux = fluxmodel(E=[3e3, 2e2], t=[1, 2, 3])
print(f'flux.shape = {flux.shape}')
print(f'flux = {flux}')

flux.shape = (1, 2, 3)
flux = [[[1.11111111e-09 1.11111111e-09 1.11111111e-09]
  [2.50000000e-07 2.50000000e-07 2.50000000e-07]]]