Katu: a fast open-source full lepto-hadronic kinetic code suitable for Bayesian inference modeling of blazars

Taking the derivative as an operator, the solution of (LABEL:kinetic_equation_compact) is:

where $\mathcal{Q}$, $A^{\prime}$ and $B^{\prime}$ all depend on time.

To obtain the first order solution to (LABEL:Exponential_Integrator_Solution), we consider that every term varies very little between $t_{0}$ and $t_{0}+\Delta_{t}$ and can be taken as constant. Then we obtain:

where the exponential and the inverse of an operator are defined through their Taylor expansions:

Note that both $A^{\prime}$ and $B^{\prime}$ depend on $\gamma$, so the cross terms of $B^{\prime}\frac{\partial}{\partial\ln\gamma}A^{\prime}$ and $B^{\prime}\frac{\partial}{\partial\ln\gamma}\left(B^{\prime}\frac{\partial}{\partial\ln\gamma}\right)$ are not null.

This step usually involves the exponentiation and inversion of an operator, which is very costly computationally. But, as we have said in the previous section, this solution is very interesting in the case that $B^{\prime}=0$ as it reduces to:

where the exponential and the inverse are easily calculable. This is the solution that we employ for neutral particles (neutrons, neutral pions and neutrinos).

An alternative way of solving (LABEL:kinetic_equation) is to discretize the derivative in energy as the first step and then look for a solution. In order to make the solution stable for high values of $\Delta_{t}$, we work from an implicit expression for the derivative. In this case, it is better to work with logarithmic derivatives, as this makes the numerical solution more stable. To account for the possible flows of particles due to the acceleration and cooling, we have to look for an equilibrium point where the two fluxes converge and separate the solution in two parts at the equilibrium energy.

We note from (LABEL:energy_change) that there is an equilibrium point in energy at $\gamma_{eq}=-\frac{1}{\tau_{acc}S}$ (remember that $S<0$) and that for $\gamma<\gamma_{eq}$, we have that $\frac{\partial\gamma}{\partial t}>0$ and conversely, for $\gamma>\gamma_{eq}$, we have that $\frac{\partial\gamma}{\partial t}<0$. This gives us the energy direction and thus, the direction in which we have to take the derivatives when discretizing. This is equivalent to upwind schemes for flows of matter, although in energy space.

In this case the discretized version of (LABEL:kinetic_equation_compact_log_log), taking index $i$ as the index of the highest energy that is lower than $\gamma_{eq}$, now reads:

Discretizing now in time, and isolating the terms in $\ln n_{j}^{t+1}$ we obtain:

For convenience, we define the following constants:

The solution can now be found in terms of the Lambert W-function:

The Lambert function has the property that $\ln W\left(x\right)=\ln x-W\left(x\right)$ for $x>0$, so we can obtain the logarithm of the population as:

By imposing some boundary conditions for $\ln n_{0}^{t+1}$ and $\ln n_{n}^{t+1}$, the rest of the populations can be solved for in a straightforward manner. This is the solution that we employ for charged particles (protons, electrons and positrons, muons and charged pions).

Katu is initialized according to a configuration file with an extensive range of options (explained in Appendix B). Here, we explain how the configuration settings are translated into the model variables.

The initial electron and proton populations are set using background density $\left(\rho\right)$, ratio of protons to electrons $\left(\eta\right)$ and the shapes of their distributions over energies as set in the configuration file. The resulting initial number densities for electrons and protons respectively are therefore given by:

Here $f_{e}$ is the normalized distribution function for electrons and $N_{e}$ the total electron number density (likewise for protons).

For the external injection of particles, we follow a similar procedure, except that instead of using the density of the fluid $\rho$, we use the total particle luminosity $L$. Then, if we inject protons and electrons in a plasma volume $\left(V\right)$ we have:

where $Q_{0,e}$ and $Q_{0,p}$ are respectively the number densities of injected electrons and protons, respectively. Injecting $\eta$ protons per electron $\left(Q_{p,0}=\eta Q_{e,0}\right)$, and calculating the integrals, we have that:

where $\left\langle\gamma\right\rangle$ is the average energy of the injected particles. Solving for $Q_{e,0}$, we can obtain that the injection distributions are given by:

The external injection of photons is much simpler as only one particle kind is involved and thus we have:

where $L_{\gamma}$ is the luminosity of the external photon field and $\left\langle\varepsilon\right\rangle$ is its normalized average energy.

The value of the free escape timescale is calculated from the volume as:

where $R$ is the radius of the sphere and $h$ is the height of the disk.

Katu evolves the population of the particles according to the kinetic equation (1). As such, it would continuing stepping indefinitely if no stopping criteria are provided. In our modelling, and in the examples that we provide with the code, we have added two stopping criteria.

In the first case, we stop the modelling if the elapsed time is bigger than a $t_{max}$ set in the configuration file.

In the second case, we stop if the simulation reaches a steady state with a tolerance better than the one required in the configuration file. For this, we check the evolution of three species of particles: electrons, photons and neutrinos, according to the following formula:

where $n$ is the number density for a given type of particle.

A number of strategies have been implemented to optimize the computational efficiency of our code.

We define a dataset $\mathcal{D}$ as the combination of observations $y$ and their associated errors $\sigma$. A model $\mathcal{M}$ is defined as some function that can be applied to a set of parameters $\theta$ to obtain some model predictions $y_{\mathcal{M}}$ that can be compared with the dataset. Then, we can define the probability that a dataset has been obtained given a model and a set or parameters as:

where $\mathcal{L}$ is called the likelihood. Likewise, we can also define the a priori probability, or belief, of the parameters that we have used:

where $\pi$ is denoted the prior. For example, when comparing leptonic and lepto-hadronic models, the parameter that most determines this division is $\eta$. If we have any a priori information suggesting a certain distribution for $\eta$, we can include this information as a prior. If not, the prior is typically set to be flat. Note however that this does not mean that no prior assumptions were included: a flat prior in logarithmic space will for example have a different impact than a flat prior in linear space.

By integrating over the whole parameter space $\Theta$, we can obtain the probability of obtaining a dataset given a model:

where $\mathcal{Z}$ is the evidence of the model, or, simplifying the notation:

The evidence can be understood as a weighted average of the likelihood $\mathcal{L}$ over the whole parameter space with the prior $\pi$ as a weight function. Thus, only the regions where the product of both the likelihood and the prior is high will mainly contribute to the value of the evidence. By integrating over the whole parameter space with a normalized prior, we are implicitly dividing by the volume of the space so models with more parameters that do not improve the likelihood will obtain a lower evidence.

By applying Bayes’ theorem, we can find the probability of a model given a dataset:

where $P\left(\mathcal{M}\right)$ and $P\left(\mathcal{D}\right)$ are the a priori probability, or belief, of the model and the data respectively. Here, $P\left(\mathcal{D}\right)$ denotes the probability of measuring the dataset $\mathcal{D}$ from the underlying physics and the observation method. $P\left(\mathcal{M}\right)$ encapsulates our prior belief in the model $\mathcal{M}$ that is not already reflected in the priors of the parameters. For example, if we know that blazars could be fit from two different models with a ratio of $1$ to $4$, we can incorporate this information in this term.

Denoting two models by $a$ and $b$, we can take the ratio of their probabilities to obtain:

Assuming that both models have the same a priori probability, we are left with the ratio of evidences as a way of comparing two models, so that if $R>1$ model $a$ is more likely and vice versa, if $R<1$, model $b$ is more likely.

As with all criteria for model selection, the actual amount of (relative) confidence in a given model that is linked to a given numerical value of R remains subjective. Taking as a guideline Kass & Raftery (1995), we give the following references for the value of $R$: $1-3.2$, ‘not worth more than a bare mention’; $3.2-10$, ‘substantial’; $10-100$, ‘strong’ and $>100$, ‘decisive’.

By inverting (LABEL:evidence_expanded) and the conditioning, we can find the probability of a point in parameter space given a dataset and a model:

which we call the posterior distribution of $\theta$: $\mathcal{P}\left(\theta\right)$. Simplifying the notation, we obtain:

For the purpose of parameter estimation, it is not necessary to compute the evidence.

Assuming that the data follow a Gaussian distribution, we can define the log-likelihood⁴⁴4The reader should not confuse the Likelihood $\left(\mathcal{L}\right)$ with the log-likelihood $\left(L\right)$. associated between a dataset and the model predictions as:

The first term is a constant for every log-likelihood, as it does not depend on the model, only on the errors of the data. It can be proven that this term appears as a constant multiplicative factor upon the calculation of the evidence, so it divides out when comparing the evidence of two models.

MCMC and Nested Sampling are two different approaches to analyse, using bayesian inference, a model given its likelihood. MCMC can sample from a given function, which in this case it is the likelihood of the model, and thus, we can reconstruct the shape of the likelihood and the posterior probabilities of any of its parameters. Nested sampling works by sampling from the likelihood function in a monotonically increasingly way. By assigning a statistical weight to each of the samples, the shape of the likelihood can be reconstructed and the evidence of the model can be calculated.

Even though sampling using MCMC is simple in principle, it has some drawbacks when the shape of the sampled distribution has various maxima and/or is very different from a simple gaussian. In order to avoid partially these problems, Goodman & Weare (2010) introduce the ’Ensemble Samplers with Affine Invariance’ and Foreman-Mackey et al. (2013) creates, using ideas from that work, emcee. The MCMC steps introduced in Goodman & Weare (2010) and implemented in emcee allows for sampling efficiently from any distribution that can be reduced to a multivariate gaussian by an affine transformation, increasing greatly the distributions for which MCMC is efficient at the cost of increasing the memory requirements. In particular, whereas simple MCMC uses a ’walker’ to traverse the parameter space, emcee uses an ensemble of them.

Even though the affine invariant ensemble samplers and emcee are a great statistical tool, they still have some problems inherited from the basic MCMC analysis. First, MCMC tends to perform very badly in situations where multiple modes exist, as it can take a long time for a walker to go from one maximum to another because in its path there is a volume of likelihood that is taken with lower probability than other steps. Second, MCMC needs some period of burn-in until the sampler starts to produce samples that are independent and representative of the sampled probability distribution. Third, not every distribution is well represented by affine transformations of a multivariate gaussian, so there are cases where emcee will still perform not optimally.

Nevertheless, after some iterations, we obtain a chain of parameters for each walker from which we can derive the sought for posterior probabilities by simply plotting the appropiate histograms.

In contrast to MCMC analysis, nested sampled is aimed towards the calculation of the evidence. The basic idea is that, by creating a sequence of points with increasingly likelihood and by assigning a weight to each of them depending on their position in the sequence, the evidence can be calculated. As the value of the likelihood increases in the sequence, the probability of drawing a new point with a better value decreases. multinest avoids this problem by using a set of points as ’live points’ and using them to define ellipses from where to drawn the new points (Feroz et al., 2009). This way, the sampling is much more efficient and the sequence can grow faster.

As the set of ’live points’ define the sampled volume, it can happen that this volume can be separated into isolated sub-volumes. In these cases, each of the sub-volumes is associated with one maximum of the likelihood function, and Nested Sampling, and multinest by extension, can detect and handle these separated maxima with ease.

Even though the length of the sequence can be made as big as wanted, there is a point where the last points contribute very little to the calculation of the evidence. multinest, by default, is set to stop the calculation of the evidence when the best point of the sequence increases the value of the log evidence by less than some tolerance.

Each of the points in the sequence is assigned a weight, and thus, can be used to create the histograms of the posterior probabilities of the parameters.

The blazar Markarian 421 (Mrk 421) has been observed for decades. Mrk 421 is one of the nearest blazars to earth, and therefore one of the brightest. Mrk 421 is a BL Lacertae object, which tend to be variable in their flux levels and can typically be modelled without a need for an external source of photons.

The full dataset that we study was gathered during a year long multi-wavelength campaing during 2009, whose details are reported on Abdo et al. (2011) and we reproduce it in 4. It is important to notice that not all the experimental points follow the same trend, even when taking into account their errors. As an example of this problem, we reproduce the X-rays part of the SED in 5, where it can be clearly seen that there are two different trends with the different instruments, Swift and RXTE, where the first shows a lower flux than the second.

Similarly, observations at optical frequencies display a spectral variability that goes beyond their (very small) observational errors. To account for these systematic deviations from the base line physics captured by our model, we take a conservative approach to the sizes of the error bars of the observations (for an alternative approach when fitting complex models to broadband transient data, see Aksulu et al. (2020)). By increasing the error bars according to the prescription below, we acknowledge that our model fit is to the generic state of the blazar SED during the period of observation rather than to its instantaneous value at a given time of observation for a specific instrument. We use

where $\sigma_{i}$ is the final error that we take, $a$ is a constant, $y_{i}$ is the experimental flux and $\sigma_{i,0}$ is the experimental error. For this analysis, we have chosen to separate the value of $a$ by band and have set $a_{optical}=a_{\gamma-rays}=0.1$ and $a_{x-rays}=0.2$. In 5 we show with bands the resulting error bands that we obtain and see that a common trend can be more easily fitted between the observations.

For the case of X-rays, the existence of two trends in the data had an important impact upon the fitting procedure and its results. The first consequence is that it creates two maxima in the likelihood function where each one is linked to each of the trends. The existence of more than one maximum can cause problems to fitting algorithms. For emcee, as we have discussed in 6.2, because it needs to have walkers in both maxima as it would be complicated for the walkers to go from one maxima to the other, and because in the degenerate case of having all walkers in one of the maxima the fitting procedure would need a very high time to ’discover’ the other maxima. For multinest the multiples maxima are not such a serious problem as it has the capability of discovering, separating and analysing each of the maxima separately, although in some cases, depending on the shape of the likelihood function this is not easy, or even doable. In both cases, what we would find is that the running time and the acceptance ratio of the fitting algorithms would be very high and low respectively.

The flux detected by SMA and most of the radio observatories should be handled carefully as they can not resolve the jet and instead detect the integrated flux from the galaxy. To take this into account, we have chosen to use the SMA point as an upper limit in our fits. Upper limits can be treated in different ways in a model fit. In this case, we have chosen to discard all model-generated SEDs that led to flux levels above the SMA value. Because the SMA upper limit in practice then automatically forces the SED to go below the radio upper limits at longer wavelengths, we can safely omit these from our fitting process without changing the outcome (confirmation of this is provided in 7.3).

As we have discussed in 4 and can be seen in B, Katu supports a high number of configuration options. To focus on one particular set of models, we choose to set the plasma volume to a sphere and to allow for all the particles to evolve. Additionally, all the initial population distributions have ben set to simple power laws.

The background density $\left(\rho\right)$ and its ratio of protons to electrons has been set to $10^{-26}g\cdot cm^{-3}$ and $1$ respectively. Note that these values have no reflection on the final steady state, as the initial distribution will be overwhelmed by the injected one.

For the electrons, we have chosen a lower cut-off for the energy range of $2$ and higher cut-off of $\max\left(\gamma_{e,max},\frac{m_{p}}{m_{e}}\gamma_{p,max}\right)$, where $\gamma_{e,max}$ and $\gamma_{p,max}$ are the values of the upper cut-offs of the electron and proton injection distributions respectively. These limits allow that electrons can cool down to very low energies and that all the electrons that are produced through photo-meson and the Bethe-Heitler processes can be injected in the distribution.

Photons are also produced by many processes, so their energy limits have to be chosen to reflect this. We have taken a lower cut-off of $10^{-12}$ and a higher cutoff of $\gamma_{e,0,max}$, which is the value of the higher energy cut-off for the electrons. We have also chosen to not inject photons from the exterior of the volume, so we set the value of the photon luminosity $\left(L_{\gamma}\right)$ as $0$.

For protons, we have taken that both, the initial and the injection distribution have the same shape, with a lower cut-off value of $\gamma_{p,min}=10$, a higher cut-off value of $\gamma_{p,max}=10^{6}$ and a power law index of $2$.

We assume that the charged particles scape slower than neutral particles, so we have set the charged to free escape ratio $\left(\sigma\right)$ as $10$.

Even though we do not know of any process that would preferentially accelerate one kind of particle over the other, we have chosen to set the ratio of protons to electrons $\left(\eta\right)$ as a free variable to further explore the hadronic capabilities of our model. In 7.4 we fix this ratio to $10^{-10}$ in order to have virtually no protons and thus, be able to simulate a fully leptonic model.

Numerically, we have chosen to set the tolerance to $10^{-6}$ and have set the number of bins for the primary particles as: $n_{e}=n_{p}=160$ and $n_{\gamma}=128$. Experience shows that these settings provide a good balance between numerical accuracy and reasonable convergence times for the simulations.

Finally, we are left with $8$ free variables, which are gathered in 2, along with the respective search ranges. $\gamma_{e,inj,min}$, $\gamma_{e,inj,max}$ and $p$ correspond to the parameters of the injected population of electrons. $B$ $\left(\text{G}\right)$ and $\delta$ are the magnetic field and the Doppler boosting parameter of the volume, respectively; $R$ $\left(\text{cm}\right)$ is the radius of the sphere we simulate; and $L$ $\left(\text{erg}\cdot\text{s}^{-1}\right)$ and $\eta$ correspond to the total luminosity of the injected particles (in the plasma frame) and its ratio of protons to electrons. As cosmological parameters, we have chosen those of Collaboration et al. (2019): $H_{0}=67.66$, $\Omega_{m}=0.3111$ and $\Omega_{\Lambda}=0.6889$.

The results of the two statistical packages that we have tested are presented in 3 and 4. The corresponding corner plots for both packages can be seen in LABEL:fig:[s]Corner-Plot-multinestand 7 respectively.

It can be clearly seen from 3 that both approaches give very similar results for the best parameters, the median and the $16$ and $84$ percentiles, which show the consistency of the analysis of both packages. Likewise, from both corner plots, LABEL:fig:[s]Corner-Plot-multinestand 7, we see that the posterior shape and the 2D correlations are consistent. Furthermore, we see in 4 that both SEDs are very close together, and indeed, their log likelihoods are $-39.78$ and $-39.33$ for multinest and emcee, respectively. For multinest, we can also quote the value of the log evidence, which is not available for emcee, $\ln\left(\mathcal{Z}\right)=-65.750\pm 0.232$.

The choice that we made in 7.1, namely to choose the SMA datapoint as an upper limit and ignore the rest of the radio detections on the basis that the SMA was the most restricting one, is correct, as none of the best fits of the SED go above any of the radio datapoints.

In computational terms, it is complicated to compare both algorithms due to how they handle points with null likelihood, as multinest silently discards them while emcee treats them as any other point; and how the first steps of emcee had been handled. multinest has taken about a day of computations and at least on the order of $120\,000$ evaluations of the log likelihood function, whereas emcee has taken more than $400\,000$ evaluations and more than double the time, although part of these evaluations were done in initial short runs and due to leaving emcee to work overnight without a clear stopping condition.⁵⁵5All the computations were done on a server with 32 available cores with multinest and emcee calculating 8 log likelihoods in parallel and Katu with 4 threads. Even though emcee seems worse from the point of view of finishing the runs, it tends to give a result for the best log likelihood much sooner, and usually better, than multinest.

From 4, it can be seen that in the range of $\sim 10^{4}-10^{7}$ eV, there is a small hump between the synchrotron and inverse Compton peaks. In our simulations, this is caused by the pairs created by the Bethe-Heitler process Mastichiadis et al. (2005); Petropoulou & Mastichiadis (2015), highlighting the need for protons in order to completely fit the SED.

Abdo et al. (2011) estimate the luminosity of the blazar between $2.6\times 10^{46}$ and $1.2\times 10^{47}$. Taking that the luminosity from the jet frame to the blazar frame transforms as $L_{blazar}\approx\delta^{2}L_{jet}/2$, we obtain $L_{blazar,multinest}\approx 6.44\times 10^{48}$ and $L_{blazar,emcee}\approx 9.98\times 10^{48}$. Both values are clearly above the estimation of Abdo et al. (2011), an issue shared with other lepto hadronic models Liodakis & Petropoulou (2020).

To verify that the small hump that appears between the peaks in 4 is indeed caused by the presence of protons and to compare our results against a purely leptonic model, we choose to repeat the statistical analysis of Mrk 421’s data with multinest but setting $\eta=10^{-10}$. This ratio of injected protons to electrons means that, effectively, the hadronic processes have no impact upon the final fit. The choice of multinest over emcee is to allow us to obtain a value for the evidence that we can compare against the evidence for the hadronic model.

As we can see in 4, the best SED that the leptonic model can produce is unable to fit the high energy part of the X-ray observations, and this is indeed reflected in the best log likelihood that it can achieve: $-42.89$. We also find that the log evidence has a value of $\ln\left(\mathcal{Z}\right)=-67.685\pm 0.231$.

Assuming that both models have the same a priori probability, we find that the hadronic model is favoured as:

The corresponding Bayes’ Factor is $6.92_{-1.93}^{+2.68}$. This value is consistent with the qualification substantial for the evidence in favour of the lepto-hadronic model, as discussed in 6.1. Note that this measure already accounts for the different number of parameters between the two models.

Lepto-hadronic models can successfully explain some of the features of blazar SEDs, such as the gamma ray peak being caused by neutral pion decay or proton synchrotron Mannheim (1993); Cerruti et al. (2015). But, leptonic models can also fit this peak with inverse compton scattering of synchrotron photons, and the energy requirements for the proton population are usually above the Eddington luminosity Liodakis & Petropoulou (2020). Some characteristics of the SEDs can only be explained by a hadronic population, such as a population of pairs produced by the Bethe-Heitler pair production process, or the presence of neutrinos Mastichiadis et al. (2005); Petropoulou & Mastichiadis (2015); Petropoulou et al. (2015); Rodrigues (2019); Rodrigues et al. (2020). However, to our knowledge, this is the first time that the dichotomy of leptonic vs. lepto-hadronic modelling has been explored from a fully bayesian point of view using kinetic equation modelling.