SOPRANO, our numerical solution of the kinetic equations, is a Python/C based numerical tool, which takes as an input the spectral injection rate of the particles (e.g., electrons and/or protons), and computes the time evolution of all relevant particles, including the secondaries, as they interact and cool, allowing the computation of the broadband emission spectrum at any given period. Our solution relies on two underlying assumptions : i) the space is homogeneous and ii) particle distribution functions are isotropic.
Here, we summarize all terms appearing in the kinetic equations for all particle species. We denote Q, S and C as the source, sink and cooling terms, respectively. The contribution of inverse Compton scattering is denoted by RIC for the photons and it is a cooling term for the leptons.
The following nonthermal radiative processes are involved in SOPRANO
All of the models allow the use of an arbitrary user-defined shape of the particle energy distribution.
Photons are produced by the synchrotron radiation of all charged particles and by the decay of neutral pions, π0. They are absorbed by pair production and redistributed in energy by inverse Compton scattering. We neglect the absorption of photons in the photo-pion and photo-pair processes. We did not consider synchrotron self-absorption and are planning to include it in the next version. The resulting kinetic equation takes the form
Radiative models
SOPRANO allows to inject (initialize) all particle flavors (e.g. proton, pions, muons, electrons/positrons, photons etc.), in either permanent or episodic regimes, and observe the evolution of those and secondary particle spectra. In its current version, SOPRANO evolves the distribution functions of the following particles:
|
The following nonthermal radiative processes are involved in SOPRANO
|
Photons are produced by the synchrotron radiation of all charged particles and by the decay of neutral pions, π0. They are absorbed by pair production and redistributed in energy by inverse Compton scattering. We neglect the absorption of photons in the photo-pion and photo-pair processes. We did not consider synchrotron self-absorption and are planning to include it in the next version. The resulting kinetic equation takes the form
where the last sum runs on all charged particles.
Leptons (electrons and positrons) are considered as a single species. They are created by muon decay, Bethe-Heitler photo-pair production and two photon recombination. They also undergo synchrotron cooling such that the final kinetic equation reads as
Protons are losing energy by synchrotron emission, photo-pion and photo-pair interactions .Protons are produced in photo-hadronic interactions between photons and neutrons, and are turned to neutrons for a substantial fraction of photo-pion interactions. The kinetic equation takes the form
Neutrons are produced in photo-pion interactions and turned to protons by the same process. The kinetic equation takes the form
In the current version of the code, we do not include neutron decay. Indeed, for the very large particle Lorentz factor involved, neutrons would escape the source before decaying. In principle, neutrinos produced by neutron decay should contribute to the observed overall signal. But since we are considering models in which the neutron number is always much smaller than the proton number, we can safely neglect this contribution.
Charged pions π+ and π−, are produced by photo-pion interactions. Then, they cool via synchrotron emission and decay. The kinetic equation for both species takes the form
The kinetic equations were solved independently for π+ and π− since the branching ratio in photo-pion production is different for negative and positive pions. This impacts the production ratio between the different neutrino species further.
Neutral pions have a kinetic equation similar to that of charged pions but without synchrotron cooling.
Muons are produced from the decay of charged pions. They lose energy by synchrotron radiation and decay. Therefore, the kinetic equation is
Muon and electron neutrinos and anti-neutrinos are produced in the decay of pions and muons. We consider the two flavours independently, but neutrino and anti-neutrinos of the same flavour are combined.
For each of the processes, the details of the terms (Q, S, C and R) together with the cross-sections used in SOPRANO are given in the paper [Gasparyan et al. 2021].
Numerical scheme
In order to numerically solve the kinetic equations presented in Radiative models, the partial differential equations are discretized in time and in energy space. There are several well-known techniques to express derivatives, dealing with diffusion types of equations. In SOPRANO, a flavour of finite volume approximation, based on the discontinuous Galerkin method, is applied. This method fits nicely for solving integro-differential equations. The prescription of finite volume allows us to conserve particle number to machine accuracy for all processes which conserve particle number (e.g. Compton scattering). Energy conservation is also enforced by specific choices for the fluxes for diffusion-like terms or redistribution of particles between adjacent energy cells. The difficulty in our numerical implementation is in the computation of the 3- to 5-dimensional integrals which approximate the rates on each energy bin can be computed once and tabulated, which improves the computational speed drastically.For the energy discretization, a numerical grid for the energy of all particles is introduced. The grid elements are equally spaced in logspace. The number of grid intervals are not the same for different particle species and are listed in Table. On each energy cell I, we approximate the distribution function by a polynomial. In the current version we restrict to first order and consider the Legendre polynomial as the basis function. Therefore, on each energy cell I, the distribution function is approximated by
where the first order Legendre polynomial on the energy cell I is
Here
are the boundaries of energy cell I and where we introduced the additional notation 
.
We seek the weak formulation of all kinetic equations on all energy intervals I. For this, we multiply both sides of the equation by L0 and integrate over I. After simplification, we obtain a set of coupled ordinary differential equations for the time evolution of the coefficients NIi, 0 .
This specific discretization and the structure of the kinetic equation allows us to retrieve a numerical method which conserves energy as well as the number of particles when required.
Table. Characteristics of the numerical grids used by SOPRANO. The cells are equally space in logarithmic scale.
Time discretization is achieved via implicit (backward) first order Euler method. Specifically, for all leptonic processes, i.e. synchrotron, IC and pair-production, the equations are treated fully implicitly, and the non-linearity of the coupled kinetic equations are solved with the Newton-Raphson method. The hadronic processes, such as the photo-pion and photo-pair interactions, a semi-implicit method is adopted, that is to say that for the evaluation of photo-pion and photo-pair collisional terms, the photon spectrum is assumed to be explicit, while the proton and neutron distribution functions are solved for implicitly. In practice, semi-implicit treatment reduces the computational cost, instead the rate of photo-pair and photo-pion interactions might be underestimated, unless the time step is carefully chosen.